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ABSTRACT 


This  thesis  presents  the  findings  of  a  study  involving  the 
pulse  testing  of  an  on-line,  double  pipe,  countercurrent,  liquid- 
liquid  heat  exchanger.  The  study  is  presented  in  three  parts;  the 
first  is  concerned  with  the  modeling  and  simulation  of  the 
exchanger,  the  second  with  the  theoretical  aspects  of  pulse 
testing,  and  the  third  with  an  experimental  pulse  testing 
investigation. 

In  the  simulation  study,  a  four-equation  model  of  the  heat 
exchanger  is  solved  as  a  lumped  parameter  system  using  a  digital 
simulation  language  (1130  CSMP) .  It  is  tested  on  the  basis  of  its 
ability  to  predict  not  only  steady  state  but  also  dynamic  conditions 
for  a  number  of  different  heat  exchangers.  The  predicted  results 
are  within  the  limits  of  accuracy  imposed  by  the  sensing  devices 
and  heat  transfer  correlations. 

As  a  prelude  to  an  experimental  investigation,  a  theoretical 
treatise  on  pulse  testing  is  given,  starting  with  a  literature  and 
mathematical  background  review.  The  latter,  in  essence,  describes 
how  frequency  response  information  may  be  obtained  by  the  Fourier 
transformation  of  time  data  from  a  pulse  test.  Three  possible 
methods  of  performing  the  Fourier  analysis  (trapezoidal,  Filon’s, 
and  fast  Fourier  transform)  are  compared  in  a  study  employing 
theoretically  generated  time  data. 
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On  the  basis  of  investigations  carried  out  by  previous  workers. 


consideration  is  given  to  the  inaccuracies  which  may  occur  in  frequency 
response  information  obtained  by  pulse  testing.  Consideration  is  also 
given  as  to  the  requirements  for  a  "good"  pulse  shape.  This  study 
culminates  in  a  set  of  guidelines  outlining  a  procedure  for  designing 
a  pulse  test. 

A  heat  exchanger  under  direct  digital  control  is  used  in  the 
experimental  work.  Computer  programs  necessary  to  implement  an  on¬ 
line  pulse,  and  to  analyze  the  resulting  time  data,  are  given.  Five 
different  pulse  shapes  are  used  and  are  shown  to  yield  amplitude 
ratio  and  phase  information  sufficiently  accurate  for  process 
characterization  work.  The  frequency  response  results  were  compared 
with  those  given  by  a  first-order  system  plus  delay,  and  by  an  over¬ 
damped  second-order  system  plus  delay.  Both  models  gave  a  sufficient 
approximation  to  the  experimental  curve  for  engineering  design  pur¬ 
poses.  The  time  constants  of  the  second-order  model  were  calculated 
on  the  basis  of  fluid  residence  times  in  the  heat  exchanger.  The 
predictive  ability  of  this  model  is  therefore  particularly  edifying 
as  it  indicates  a  very  convenient  engineering  approach  for  describing 
the  dynamics  of  heat  exchangers. 
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1  INTRODUCTION 


The  first  part  of  this  thesis  is  concerned  with  the  modeling 
and  simulation  of  a  double  pipe  ,  countercurrent ,  liquid-liquid  heat 
exchanger . 

Prior  to  the  development  of  a  model,  it  was  necessary  to  deter¬ 
mine  satisfactory  heat  transfer  correlations  for  the  tube,  tube-to- 
annulus,  and  shell  film  coefficients.  A  review  of  the  literature 
indicated  that,  for  the  tube  film  coefficient  and  shell  film  coeffici¬ 
ent,  the  Dittus-Boelter  equation  and  its  modified  form  were  recommended 
by  most  authors.  However,  for  the  calculation  of  the  tube-to-annulus 
film  coefficient,  correlations  by  Stein  and  Begell  (53)  and  by  McAdams 
(40)  have  been  suggested.  Subsequent  simulation  work  indicated  that  the 
McAdams'  correlation  gave  slightly  better  results  and  was  therefore 
used  in  later  studies. 

The  mathematical  model  for  the  heat  exchanger  was  obtained  by 
deriving  the  unsteady  state  heat  balances  for  the  two  fluids  and  two 
tubes.  The  four  nonlinear  partial  differential  equations  which 
resulted  were  used  as  the  basis  of  the  simulation  work  and  are  referred 
to  as  the  general  four-equation  model. 

The  general  four-equation  model  defies  a  straightforward  analytic 
solution,  and  a  literature  survey  showed  that  previous  workers  have 
resorted  to  a  number  of  simplifying  assumptions  before  attempting  a 
solution.  Having  simplified  the  model,  various  techniques  have  been 
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used  to  solve  the  equations  and  thus  determine  the  transient  solution. 

The  three  basic  techniques  used  are  Laplace  transformation,  lumped 
parameter,  and  finite  difference  approximation  methods. 

In  recent  years  the  development  of  digital  simulation  languages 
has  provided  a  convenient  method  of  solving  lumped  parameter  problems. 
Advantage  was  taken  of  this  technique  and  the  four  equations  used  in 
this  study  were  solved  as  a  four-lump  model,  using  the  Continuous 
Systems  Modeling  Program  (CSMP)  available  with  the  Department  of 
Chemical  and  Petroleum  Engineering's  IBM  1800  computer.  One  major 
improvement  over  previous  methods  of  solution  was  that  the  nonlinear 
nature  of  the  model  was  preserved.  The  film  coefficients,  primary  cause 
of  the  equations  being  nonlinear,  were  solved  by  special  subroutines 
written  by  the  author  and  incorporated  into  the  CSMP  program. 

The  general  four-equation  four-lump  model  was  thoroughly  tested, 
the  basis  of  the  testing  being  its  ability  to  predict  the  steady  state 
and  transient  temperature  values  for  a  number  of  given  heat  exchanger 
systems.  Fisher  (23)  presents  the  results  of  several  steady  state  runs 
using  various  heat  exchanger  configurations.  The  model  was  initially 
tested  by  simulating  some  of  his  runs.  Further  steady  state  simulation 
work  was  carried  out  using  the  same  double  pipe,  countercurrent,  liquid- 
liquid  heat  exchanger  as  used  by  Lees  (35) .  As  a  final  test  of  the 
model,  three  runs  were  made  in  which  the  test  heat  exchanger  was  subjected 
to  various  tube  flow  disturbances.  These  runs  were  simulated  and  a  com¬ 
parison  made  of  the  actual  temperature  transients  versus  those  predicted 
by  the  model.  The  ability  of  the  model  to  simulate  the  steady  state  as 
well  as  transient  runs  was  found  to  be  within  the  accuracy  limitations 
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imposed  by  the  heat  transfer  correlations.  This  indicates  the  potential 
of  solving  nonlinear  partial  differential  equations  from  a  lumped 
parameter  approach  by  digital  simulation. 
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2  HEAT  TRANSFER  CONSIDERATIONS 

2.1  Heat  Transfer  Correlations 

To  determine,  for  a  double  pipe  heat  exchanger,  the  rate  of  heat 
transfer  from  either  the  tube  or  shell  to  the  annular  fluid,  the  pertin¬ 
ent  heat  transfer  film  coefficients  must  be  known.  These  coefficients 
are  affected  by  so  many  system  parameters  that  it  has  been  found 
expedient  to  arrange  these  terms  in  a  series  of  dimensionless  groups  (14). 

By  applying  the  principles  of  dimensional  analysis  to  the 
problem  of  heat  transfer  by  convection,  the  following  general  relation¬ 
ship  is  obtained: 

Nu  =  f  [  (Re)  (Pr)  (Gr)  ]  ...  2-1 

where  Nu  =  h  d  =  Nusselt  number 

~ 

Re  =  v  d  p  =  Reynolds  number 

y 

Pr  =  Cp  y  =  Prandtl  number 

k 

Gr  =  B  g  AT  d3  p2  =  Grashof  number 

u2 

In  the  case  of  natural  convection,  the  fluid  velocity  is  solely 
a  function  of  the  buoyancy  effects.  The  variation  of  density  with 
temperature,  as  represented  by  the  coefficient  of  thermal  expansion 
(B),  is  expressed  in  the  Grashof  number.  The  importance  of  this  group 
is  such  that  the  Reynolds  number  may  be  omitted. 
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For  forced  convection,  the  fluid  velocity  is  high  enough  that 
fluid  motion  due  to  buoyancy  effects  is  negligible,  and  the  Reynolds 
number  is  used  to  the  exclusion  of  the  Grashof  number.  Thus: 

for  natural  convection  Nu  =  f  [  (Pr)  (Gr)  ] 
for  forced  convection  Nu  =  f  [  (Re)  (Pr)  ] 

In  this  work,  only  turbulent  systems  are  to  be  studied  and  heat 
transfer  by  natural  convection  will  not  be  discussed  further. 


Probably  the  best  known  empirical  correlation  for  forced  convec¬ 
tion  heat  transfer  in  tubes  is  that  presented  by  Dittus  and  Boelter  (15) : 

(Nu),  =  0.023  (Re)°‘8  (Pr)"  . . .  2- 

b  b  b 

where  subscript  b  indicates  that  the  fluid  properties  are  evaluated  at 
the  arithmetic  mean  bulk  temperature.  For  this  correlation,  the 
following  conditions  apply: 


(i)  Re  No.  >  10,000 

(ii)  Viscosity  not  more  than  twice  that  of  water 

(iii)  n  =  0.4  for  heating,  n  =  0.3  for  cooling 

(iv)  L,  >  60 
d 


Equation  2-2  is  not  applicable  for  fluids  with  high  Prandtl  numbers. 
Seider  and  Tate  (50)  proposed  an  equation  which  is  applicable  to  those 
fluids  which  experience  large  viscosity  changes  for  small  temperature 
changes : 


(Nu)b  =  0.027^u 


t) 


0.14 


(Re) 


0.8 


(Pr) 


0.33 


. ..  2-3 


where  p  is  the  viscosity  at  the  mean  bulk  temperature  and  Ps  the 
viscosity  at  the  mean  surface  temperature  of  the  wall.  This  work 


' 
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employs  only  water  as  the  heat  transfer  medium  and  thus  correlations 
of  the  form  given  in  Equation  2-3  will  not  be  used. 


For  the  case  of  a  double  pipe  heat  exchanger,  a  difficulty 
arises  in  the  representation  of  the  two  diameters  in  the  film  coeffici¬ 
ent  correlations.  Monrad  and  Pelton  (41),  whose  investigation  covers 
heat  transfer  for  diameter  ratios  of  1.65,  2.45  and  17,  used  water  and 
air  over  a  temperature  range  of  80  degrees  F  to  180  degrees  F  and 
Reynolds  numbers  from  12,000  to  220,000.  They  found  that  a  modified 
Dittus-Boelter  equation  satisfactorily  represented  (within  ±  10  per 
cent)  the  case  of  heat  transfer  from  shell  to  annulus: 


and  d 
d 


2 

T 


0.023 


4  x  (hydraulic  radius) 


inside  diameter  of  the  shell 
outside  diameter  of  the  tube. 


. ..  2-4 


For  heat  transfer  from  tube  to  annulus,  Monrad  and  Pelton 
suggested : 


In  1945,  Wiegand  (58)  reviewed  the  literature  up  to  that  time  and 
endorsed  Monrad  and  Pelton ’s  correlation  for  the  shell-to-annulus 
coefficient.  However,  for  the  tube-to-annulus  coefficient  he  proposed 
that  their  (d2/dT)  factor  be  raised  to  the  0.45  power,  to  agree  with 
the  work  of  a  number  of  other  authors. 


.(  •  \OSS  o3  000 « SI  raoTfl  aisdniun  sblonystf 
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In  1958  Stein  and  Begell  (53)  carried  out  an  intensive  study  of 
the  problem  of  forced  convective  heat  transfer  in  annuli.  Nearly  900 
values  of  local  heat  transfer  coefficients  were  correlated  for  water 
flowing  through  long  annuli  1/8-inch,  1/4-inch  and  3/8-inch  wide.  All 
heat  transfer  coefficients  were  computed  for  (L/d^)  ratios  larger  than 
150.  They  found  the  most  accurate  correlation  (±  10  per  cent)  for  the 
tube-to-annulus  coefficient  to  be: 


(Nu)f  = 

where  t^ 

and  t, 
b 

t 

s 


0.020  de  ^Cp  \i 


=  average  bulk  temperature  of  the  fluid 

=  average  surface  temperature  on  the  outside  of  the 
tube  wall. 


McAdams  (40),  however,  recommends  the  Dittus-Boelter  equation, 
modified  to  take  into  account  the  equivalent  diameter,  i.e.: 


where  n  =  0.4  for  heating,  n  =  0.3  for  cooling,  and  subscript  b  indi¬ 
cates  that  the  fluid  properties  are  to  be  evaluated  at  the  bulk 
temperature  of  the  fluid. 

Simulation  of  the  experimental  data  on  double  pipe  heat 
exchangers  reported  by  Fisher  (23)  showed  that  Equation  2-7  came  the 
closest  to  predicting  the  experimental  steady  state  fluid  outlet 
temperatures.  The  modified  Dittus-Boelter  equation  was  therefore 


used  in  the  remainder  of  the  study. 


•  Ybuia  aria  lo  i»bn  Lr/h^-*  s/13  ni  b&au 


In  summary,  Equation  2-2  was  used  for  the  tube  film  coefficient. 
Equation  2-7  for  the  tube-to-annulus  film  coefficient,  and  Equation  2-4 
for  the  shell  film  coefficient. 


2.2  Calculation  of  Film  Coefficients 


To  establish  the  magnitude  of  the  film  coefficients  to  be 
expected  in  this  investigation,  values  were  computed  for  a  velocity 
range  of  1.0  to  5.0  ft /sec  and  a  temperature  range  of  50  degree  F  to 
200  degrees  F.  The  correlations  used  were  those  indicated  in  the 
previous  section,  namely 

(i)  Tube  film  coefficient: 

(Nu)b  =  0.023  (Re)°'8  (Pr)£  ...  2-8 


(ii)  Tube-to-annulus  film  coefficient: 


(Nu)fe 


(iii)  Shell  film  coefficient: 


. ..  2-10 


In  this  study,  cold  water  was  being  heated  by  hot  water  in  the  annulus, 
thus,  for  Equation  2-8  n  =  0.4,  and  for  Equations  2-9  and  2-10  n  =  0.3. 


.  rmaioiitl'aoo  ilsiia  sri:J  ?Q/i 
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To  use  the  above  correlations  it  is  necessary  to  evaluate  the 
physical  properties  of  the  fluid  at  some  specified  film  or  bulk 
temperature.  The  change  in  specific  heat  and  density  of  water  over 
the  temperature  range  of  interest  is  negligible;  however,  this  is  not 
the  case  for  viscosity  and  thermal  conductivity.  A  convenient 
expression  for  the  change  in  viscosity  with  temperature  is  given  by 
Bingham  (1) : 

1  =  2.1482  {  (T  -  8.435)  +  [8078.4  +  (T  -  8.435)2]^  }  -  120  ...  2-11 

U 

where  T  =  degrees  C 

y  =  poises. 

A  suitable  expression  for  thermal  conductivity  was  not  found, 
although  numerous  tables  of  data  were  located.  A  Lagrangian  inter¬ 
polating  polynomial  was  fitted  to  the  data  given  by  Kreith  (34)  and  the 
following  equation  resulted: 

k  =  (0.00266  T3  -  3.2  T2  +  1073.33  T  +  286,000  )  x  10"6  ...  2-12 

where  T  =  degrees  F 

k  =  Btu/hr  ft  degrees  F. 

Using  the  above  expressions  the  three  film  coefficients  were 
calculated  for  water  for  two  sets  of  conditions;  the  first,  a  fluid 
velocity  of  1  ft/sec  and  50  degrees  F,  the  second,  a  fluid  velocity  of 
5  ft/sec  and  200  degrees  F.  These  two  sets  of  conditions  represented 
the  range  of  interest  of  fluid  flow  and  temperature  that  could  be 
expected  in  the  simulation  work. 


, 
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The  degree  of  variation  of  the  three  film  coefficients,  shown  in 
Table  1,  indicated  that  any  model  developed  for  the  simulation  studies 
should  incorporate,  if  possible,  the  full  correlation  for  each  film 
coefficient . 


TABLE  1 

Calculation  of  Heat  Transfer  Film  Coefficients 


Heat  transfer 
film  coefficients 

h  x 

(Btu/hr 

10  3 
ft2  °F) 

v  =  1  ft/sec 

T  =  50°  F 

v  =  5  ft/sec 

T  =  200°  F 

Tube  film 

0.245 

1.76 

Tube-to-annulus  film 

0.237 

2.01 

Shell  film 

0.202 

1.71 

' 
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3  DERIVATION  OF  MATHEMATICAL  MODEL 

The  mathematical  model  for  a  concentric  pipe  heat  exchanger  is 
obtained  by  deriving  the  unsteady  state  heat  balances  for  the  two  fluids 
and  two  pipes.  A  section  of  the  exchanger  is  shown  below: 


Fig.  1  Section  of  Heat  Exchanger 


The  basic  assumptions  employed  in  this  derivation  are: 

(i)  The  heat  exchanger  is  well  insulated,  and  there  is  no  heat 
loss  from  the  shell  to  the  environment. 

(ii)  Fluid  and  metal  densities,  and  specific  heats,  are  constant. 

(iii)  Fluid  velocity  and  temperature  profiles  are  constant 
across  any  cross-section  normal  to  the  direction  of  flow, 
and  the  radial  component  of  velocity  is  zero. 


- 
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(iv)  Axial  heat  conduction  in  both  walls  and  fluids  is 
negligible  compared  to  heat  transported  by  bulk  flow. 

(v)  The  thermal  conductivity  of  the  metal  walls  is  high 
enough  that  radial  temperature  gradients  in  the  walls 
are  negligible. 

Heat  balance  on  tube  fluid: 


VlAlplCplTl 


x 


-  v1a1p1cp1t1 


x  +  Ax 


+  h1P1Ax  (Tt  -  T1) 


3  (A1Axp1Cp1T1) 
3  t 


...  3-1 


Division  of  Equation  3-1  by  A^Axp^Cp^  gives 


VlAlplCplTl 


x 


-  viAiPiCPiTi 


x  +  Ax 


Vi  (TT  -  V 


+ 


A^Axp^Cp. 


A1p1Cp. 


=  3Ti 


...  3-2 


3t 


Taking  the  limit  of  Equation  3-2  as  Ax  ■>  0  gives 


8T1  =  -  V18T1  +  hipi  (tt  ti 

3t  3x  A^p^Cp^ 


. ..  3-3 


y 


6 


■ 
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Heat  balance  on  tube  wall: 


hTPTAx  (T2  -  Tt)  -  h1P1Ax  (Tt  -  Tx) 


|r  (V^T^tV 

o  t 


3-4 


Division  of  Equation  3-4  by  A^Axp^Cp^  gives 

3TT  =  Vt  (T2  "  V  -  hlPl  (TT  “  V 

3t  ATpTCpT  ATpTCpT 


. ..  3-5 


Heat  balance  on  annular  fluid: 


v2A2p2Cp2T2 


x  +  Ax 


-  v2A2P2Cp2T2 


X 


-  hTPTAx  (T2  -  Tt) 


h2P2Ax  (T,  -  Ts) 


9  (A  Axp  Cp  T  ) 

at  z 


...  3-6 


Division  of  Equation  3-6  by  A2Axp2Cp2  and  simplifying  gives 


3T2  =  V2  9T2  -  h2P2  (T2  “  V  -  hTPT  ^2  "  V 

9t  8x  A2P2Cp2  ^2P2^^2 


. ..  3-7 


Note  that  the  velocity  terms  in  Equations  3-3  and  3-7  are  of  opposite  sign 


since  this  is  a  countercurrent  flow  model. 
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Heat  balance  on  shell  wall: 

Assuming  negligible  heat  loss  to  the  surroundings. 


h2P2Ax  (T2  -  T  ) 


9_ 

3t 


(-AsAxpsCpsTS') 


.  ..  3-8 


Division  of  Equation  3-8  by  A  Axp  Cp  gives 

u  u  u 


3Tg  .  h2P2  (T2  -  Ts) 

AspsCps 


. ..  3-9 


Equations  3-3,  3-5,  3-7  and  3-9  constitute  the  general  four-equation 


model  studied  in  this  thesis. 
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4  REVIEW  OF  HEAT  EXCHANGER  MODEL 


In  the  previous  section,  the  four  differential  equations  which 
constitute  the  model  of  a  double  pipe,  countercurrent,  liquid-liquid 
heat  exchanger  were  derived.  These  equations  are: 


-  -  Ii!Zi  +  hipi  (tt  -  V 

it  3x  AlplCPl 


=  Vi  (T2  -  V  -  hlPl  <TT  -  V 

ATpTCpT  atPtcPt 


^2  =  v2  ^2  -  h2P2  (T2  -  V 

9t  dx  A2p2Cp2 


hP 
T  T 

A2P2Cp2 


) 


3Tg  .  h2P2  (T2  -  Tg) 

at  AgPgCpg 


In  the  assumptions  employed  in  the  derivations,  no  restrictions 
were  placed  on  the  form  of  the  heat  transfer  coefficients.  As  previously 
shown,  the  coefficients  are  functions  of  fluid  velocity  and  the  physical 
properties  of  the  system.  This  concept  is  commonly  expressed  as 

Nu  =  f  (Re,  Pr)  or  one  may  write 
h  =  f  (velocity,  temperature). 

The  coefficients  in  the  model,  therefore,  are  functions  of  their 
respective  stream  temperatures.  This  makes  the  equations  nonlinear  and 


' 


. 
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renders  the  model  practically  impossible  to  solve  analytically.  In 
order  to  circumvent  this  problem  and  facilitate  a  theoretical  analysis, 
previous  workers  have  resorted  to  various  simplifying  assumptions.  One 
way  to  reduce  the  analytical  difficulties  is  to  keep  fluid  velocities 
constant  and  only  perturb  the  temperatures.  This  is  the  reason  for  most 
of  the  early  work  on  exchanger  dynamics  being  on  temperature  forcing 
instead  of  flow  forcing.  Other  simplifying  assumptions  most  commonly 
used  are  as  follows: 

(i)  The  four-equation  model  may  be  reduced  to  a  three- 
equation  model  by  ignoring  the  effects  of  the  shell  wall 
dynamics . 

(ii)  Further  reduction  in  the  number  of  equations  is  attained 
by  neglecting  the  heat  capacity  of  the  tube  wall;  this 
results  in  a  two-equation  model. 

(iii)  A  one-equation  model  which  considers  only  one  component, 
usually  the  tube  fluid,  may  be  acceptable  in  some  cases. 
For  example,  condensing  vapor  in  the  shell,  or  a  very 
high  mass  flow  in  the  shell  of  a  short  exchanger,  could 
result  in  a  practically  constant  tube  wall  temperature. 

(iv)  The  dependence  of  the  heat  transfer  coefficients  (h)  on 
velocity  and  temperature  causes  nonlinearities  in  the 
equations.  This  may  be  avoided  by  assuming  the  coeffici¬ 
ents  to  be  constant  over  the  experimental  range  of 
interest,  and  combining  them  into  a  constant  overall  heat 
transfer  coefficient  (U) . 

(v)  A  more  rigorous  approach  than  (iv)  is  to  assume  that  the 
coefficients  are  dependent  either  on  velocity  or 
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temperature,  and  use  some  statistically  determined 
function  based  on  available  experimental  data. 
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(vi)  The  partial  differential  equations  may  be  reduced  to 
ordinary  differential  equations  by  assuming  that  the 
temperature  observed  in  both  shell  and  tube  is  constant 
with  respect  to  x,  over  some  distance  Ax.  This  "lumped 
parameter"  approach,  since  it  entails  the  solving  of 
ordinary  differential  equations,  is  very  amenable  to 
analog  simulation.  Obviously,  the  smaller  Ax  is,  and 
thus  the  more  lumps  used,  the  better  the  representation 
of  the  model  to  be  expected. 

Having  simplified  the  model,  various  techniques  have  been  used 
to  solve  the  equations,  and  thus  determine  the  transient  solution. 

Three  basic  techniques  have  been  used  by  previous  investigators.  These 
are  as  follows: 

(i)  Linearizing  the  simplified  equations  by  perturbation 

techniques,  then  applying  Laplace  transform  methods  to 
obtain  the  required  transfer  functions.  This  is  accom¬ 
plished  by  expressing  the  variables  as  the  sum  of  a  steady 
state  value  and  a  perturbation  from  steady  state. 
Substitution  of  these  values  into  the  equations,  and  drop¬ 
ping  the  terms  involving  the  product  of  two  perturbations 
gives  the  equations  in  terms  of  steady  state  and  dynamic 
quantities.  Subtraction  of  the  steady  state  equations 
yields  the  dynamic  equations,  which  have  constant 
coefficients.  Laplace  transformation  methods  can  now  be 
applied  to  produce  the  desired  transfer  function. 
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(ii)  Using  a  lumped  parameter  representation  of  the  model, 
analog  computers  can  be  used  to  solve  the  equations. 

(iii)  Finite  difference  approximations  may  be  substituted  for 
the  partial  derivatives  in  the  model.  This,  in  effect, 
sets  up  a  two-dimensional  grid  with  a  spatial  variable 
and  a  time  variable.  The  continuous  transient  tempera¬ 
ture  distribution  can  be  approximated  by  the  four 
temperatures  at  each  of  the  grid  points,  or  nodes.  Using 
the  initial  temperature  distribution,  which  is  assumed  to 
be  known,  and  the  appropriate  finite  difference  approxima¬ 
tion,  a  new  temperature  profile  can  be  calculated. 
Apparently,  there  are  dangers  of  nonconvergence  and 
unstable  solutions  associated  with  this  method. 

This  general  classification  of  simplifying  assumptions  and  solving 
techniques  is  meant  as  an  aid  to  the  reader,  for  easier  evaluation  of 
the  work  of  previous  investigators  who  are  mentioned  in  the  literature 
survey  which  follows.  This  investigation  is  not  concerned  with  tempera¬ 
ture  forcing,  thus  only  work  on  flow  forcing  systems  has  been  reviewed. 

Mozley  (43)  approaches  the  problem  of  representing  a  double  pipe 
heat  exchanger  entirely  from  a  lumped  parameter  basis.  He  assumes,  in 
the  derivation,  perfect  mixing  in  both  lumped  sections  of  the  tube  and 
annulus,  negligible  wall  capacitance,  and  the  use  of  a  constant  overall 
heat  transfer  coefficient.  Even  though  all  of  the  experimental  work, 
and  most  of  the  theoretical  analysis,  was  concerned  with  temperature- 
forced  cases,  Mozley  was  the  first  to  use  the  lumped  parameter  approach 


for  the  flow-forced  case. 
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Edwards  (19)  was  one  of  the  first  to  investigate  the  dynamics  of 
a  flow-forced  double  pipe  heat  exchanger.  He  obtained  experimental 
frequency  response  data  and  compared  this  with  a  theoretical  transfer 
function  derived  from  a  linearized  distributed  parameter  model.  A 
comparison  was  also  made  of  the  experimental  data  with  harmonic  analysis 
results  based  on  a  lumped  parameter  model.  It  was  concluded  that  the 
proposed  lumped  parameter  model  fitted  the  experimental  data  better 
than  the  proposed  distributed  parameter  model. 

Stermole  (54)  also  did  an  experimental  and  theoretical  study  of 
a  double  pipe  heat  exchanger.  As  a  basis  for  the  investigation 
frequency  and  transient  response  analysis  was  employed.  Three  models 
were  used  in  his  theoretical  work;  a  three-equation  model  with  constant 
coefficients,  a  two-equation  model  with  constant  coefficients,  and  a 
lumped  parameter  model  utilizing  two  equations.  It  was  concluded  that 
the  effects  of  wall  capacitance  and  heat  transfer  coefficient  variation 
on  the  normalized  frequency  response  results  were  minor.  It  was  also 
concluded  that,  for  simulation  purposes,  variable  coefficient  models 
should  be  used  to  represent  distributed  parameter  heat  exchange  systems 
when  large  flow  perturbations  are  expected  (±  25  per  cent  of  the  mean 
flow  rate) . 

Isom  (33),  in  his  study,  used  a  distributed  parameter  two- 
equation  model  with  an  overall  heat  transfer  coefficient.  The  fit  of 
the  distributed  parameter  solutions  to  the  experimental  data  was  poor 
and  it  was  suggested  that  inclusion  of  the  effect  of  tube  wall  capa¬ 
citance  should  give  an  improvement.  In  a  lumped  parameter  representa¬ 
tion,  based  again  on  only  the  two-equation  model,  the  overall  heat 
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transfer  coefficient  was  allowed  to  vary  with  flow,  according  to  an 
equation  of  the  form 
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U  =  a  4-  b  w 
c  +  d  w 

where  a,  b,  c  and  d  are  constants  and  w  is  the  forcing  fluid  flow  rate. 
This  model  provided  a  better  fit  of  the  experimental  amplitude  ratio 
and  phase  angle  curves.  Runs  were  made  with  four,  six  and  eight  lumps, 
and  the  difference  between  the  predicted  results  was  found  to  be  slight. 

Herron  (27)  presented  a  number  of  models  of  varying  complexity; 
the  least  restrictive  being  a  three-equation  tube  wall  capacitance 
representation.  The  emphasis  of  the  work  was  on  digital  simulation, 
and  in  his  "variable  coefficient  model"  the  heat  transfer  coefficient  is 
assumed  to  vary  with  the  0.8  power  of  the  shell  side  velocity.  On  the 
basis  of  a  two-equation  model  with  a  constant  overall  heat  transfer 
coefficient,  Herron  also  did  an  analog  computer  simulation,  and  con¬ 
cluded  that  digital  simulation  of  nonlinear  models  may  not  always  be 
different  enough  from  the  results  of  simpler  models  to  warrant  the 
increased  calculation  time.  It  was  suggested  that,  if  severe  forcing 
is  not  anticipated,  the  simpler  lumped  parameter  simulation  may  be 
sufficient . 

Privott  (46) ,  also  concerned  with  the  analysis  of  flow-forced 
concentric  tube  heat  exchangers,  investigated  the  errors  associated 
with  the  digital  simulation  of  nonlinear  systems.  His  nonlinear  model 
includes  the  dynamics  of  both  shell  and  tube  walls  (four  equations) , 
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and  the  heat  transfer  coefficients  are  functions  of  their  pertinent 
fluid  stream  velocities  to  the  0.8  power.  In  another  model,  the 
temperature  dependence  of  the  heat  transfer  coefficients  is  accounted 
for  by  using  the  fact  that,  over  moderate  temperature  ranges,  the 
heat  transfer  coefficients  for  water  can  be  approximated  as  linear 
functions  of  temperature.  For  a  given  stream  velocity,  the  heat 
transfer  coefficients  were  represented  by  the  form 

h  “  hR  [1  +  SR  (Ia  -  V  1 

where  the  subscript  R  refers  to  conditions  at  a  reference  temperature, 
the  subscript  a  indicates  the  arithmetic  average  between  the  inlet  and 
exit  values  of  the  particular  stream,  and  S  is  the  slope  of  the  linear 
approximation.  It  was  found  that  the  best  theoretical  linear  representa¬ 
tion  of  the  process  was  the  model  derived  from  the  nonlinear  model  by 
perturbation  analysis.  However,  this  model  predicted  temperature 
distributions  which  were  significantly  different  from  the  true  distribu¬ 
tions  for  even  modest  velocity  changes.  It  was  concluded  that  errors 
introduced  by  linearization  in  the  steady  state  temperature  distribu¬ 
tions  after  step  changes  in  velocity  are  strong  functions  of  the  system 
parameters,  as  well  as  the  velocity  changes. 

Fisher  (23) ,  in  his  work  on  frequency  response  of  flow-forced 
double  pipe  heat  exchangers,  carried  out  an  experimental  investigation 
on  various  sized  heat  exchangers.  A  valuable  amount  of  experimental 
data  is  presented.  Linearized ' transfer  functions  were  derived  for  the 
frequency  response  of  both  fluid  outlet  temperatures  to  disturbances 
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in  either  flow  rate  under  countercurrent  flow  conditions.  The 
transfer  functions  were  based  on  a  three-equation  model  and  included 
the  effects  of  tube  wall  heat  capacity  and  variations  in  heat  transfer 
coefficients  with  flow.  The  mathematical  models  accurately  predicted 
the  frequency  response  of  the  larger  exchangers  operated  at  high 
Reynolds  numbers.  It  was  noticed,  however,  that  as  the  heat  exchanger 
size  was  reduced,  the  attenuation  and  phase  shift  of  the  frequency 
response  of  the  fluid  outlet  temperatures  increased,  and  became 
significantly  greater  than  that  theoretically  predicted.  This  was 
attributed  to  axial  mixing  effects. 
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5  CSMP  SIMULATION 


5.1  Introduction 

In  recent  years  digital  simulation  languages  have  become 
available,  enabling  lumped  parameter  problems  to  be  solved  digitally. 

The  four-equation  model  used  in  this  study  was  solved  in  this  manner, 
using  the  1130  Continuous  Systems  Modeling  Program  (CSMP)  available 
on  the  IBM  1800  computer  in  the  Data  Acquisition,  Control  and 
Simulation  (DACS)  Centre  of  the  Department  of  Chemical  and  Petroleum 
Engineering . 

CSMP  (31,  32)  is  a  "digital  analog  simulator  program"  in  which 
the  functional  blocks  of  the  input  language  represent  the  elements  and 
organization  of  an  analog  computer.  It  includes  a  complement  of 
standard  simulation  elements  such  as  integrators,  multipliers  and 
summers,  as  well  as  a  group  of  special  elements  such  as  dead  space, 
limiter,  and  zero  order  hold.  Each  element  is  provided  with  a 
diagrammatic  symbol  and  a  language  symbol.  The  user  starts  by 
developing  a  block  diagram  representation  of  his  process  in  much  the 
same  way  as  an  analog  computer  representation.  The  diagram  showing 
the  elements  and  their  interconnections  is  then  translated  into  a  set 
of  CSMP  language  statements. 

CSMP  has  the  very  useful  facility  of  being  able  to  incorporate 
user-supplied  subroutines.  This  was  taken  advantage  of  by  writing  three 
special  subroutines  to  calculate  the  three  heat  transfer  coefficients, 
as  given  by  the  correlations  discussed  in  Chapter  2.  A  listing  of 
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these  subroutines,  plus  a  subroutine  to  produce  a  pulse,  is  given  in 
Appendix  I-A. 

Thus,  this  method  of  solution  not  only  used  the  four-equation 
model,  but  also  obviated  the  need  to  simplify  the  heat  transfer 
correlations,  as  previous  investigators  had  done.  A  total  of  69 
elements  was  used  to  represent  the  four-equation  model  as  a  series  of 
four  lumps.  A  diagram  representing  one  of  the  lumps  is  given  in 
Figure  2,  and  a  listing  of  the  CSMP  configuration  and  specification 
statements  is  given  in  Appendix  I-B. 

As  previously  mentioned,  two  possible  correlations  were 
available  for  the  tube-to-annulus  heat  transfer  film  coefficient, 
one  by  McAdams,  the  other  by  Stein  and  Begell.  Both  correlations  were 
used  in  the  initial  simulation  studies  and  the  results  compared. 


5.2  Simulation  Study  using  the  Experimental  Data  of  Fisher 

In  the  course  of  his  dynamic  studies,  Fisher  (23)  carried  out 
an  experimental  investigation  on  various  sized  heat  exchangers. 
Steady  state  conditions  for  four  of  his  experimental  runs  were 
simulated,  and  the  model  used  to  predict  the  two  recorded  fluid 
outlet  temperatures.  Predictions  of  these  outlet  temperatures  were 
made  with  the  CSMP  model  using  the  tube-to-annulus  correlation  as 
given  by  McAdams,  and  then  using  the  correlation  given  by  Stein  and 
Begell.  The  results  are  shown  in  Table  2. 
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FIG. 2  C5MP  MODEL-  i”  LUMP 


Legend  for  Figure  2 


Parameters 


P^  =  ^  at  steady  state 

P0  =  -  1.0 

P8 

=  -  PT 

A2P2Cp2 

Ax 

Pq 

=  -  P2 

P3  =  P1 

y 

A2P2Cp2 

axp 1cp1 

pio 

=  T  „  at  steady  state 
b  ,  1 

P4  =  i  at  steady  state 

pn 

=  P2 

*5  •  -  Pl 

AspsCps 

ATpTCpT 

P-,  o 

=  v2 

P 

P6  =  T 

12 

Ax 

atptcpt 

Py  =  ^21  at  steady  state 

P13 

=  -  V2 

Ax 

Subscripts 

1  -  Refers  to  the  tube  liquid 

T  - 

refers  to  the  tube 

2  -  Refers  to  the  shell  liquid  S  - 

refers  to  the  shell 

Block  Symbols 

I  -  Integrator 

+  - 

Summer 

X  -  Multiplier 

W  - 

Weighted  summer 

SUB^  -  Special  subroutine  for 

calculating 

the  tube  film  coefficient 

SUB,  -  Special  subroutine  for 
coefficient 

calculating 

the  tube-to-annulus  film 

SUB^  -  Special  subroutine  for 

calculating 

the  shell  film  coefficient 

'  -  .'Ol 
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Results  of  CSMP  Simulation  employing  the  Experimental  Data  reported  by  Fisher  (23) 
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An  explanation  of  the  terms  used  in  the  table  is  as  follows: 


F.D. 


S.B. 

M 

Vqt 


Experimental  data  of  Fisher 

Results  using  Stein  and  Begell  correlation 

Results  using  McAdams'  recommended  correlation 

Ratio  of  heat  transfer  between  the  two  fluid  streams; 

subscript  S  refers  to  the  shell  fluid,  subscript  T 

refers  to  the  tube  fluid. 

Cross-sectional  area  of  the  tube 
Cross-sectional  area  of  the  annulus 


The  heat  exchanger  configurations  for  each  run  are  given  below: 


Fisher’s 


Run  No. 

Heat  Exchanger  Size 

at  (ft2) 

As  (ft2 

2 

0.75"  x  1.125"  x  6’ 

0.00230 

0.00233 

20 

0.50"  x  0.75"  x  4' 

0.00100 

0.00094 

23 

0.50"  x  0.75"  x  4' 

0.00100 

0.00094 

25 

0.375"  x  0.625"  x  3' 

0.00597 

0.00743 

Simulation  of  Runs  2  and  25  gave  the  best  prediction  of  tube  and 
shell  outlet  temperatures.  In  Run  2,  both  heat  transfer  correlations 
gave  predicted  outlet  shell  temperatures  to  within  0.5  degrees  F  and 
both  tube  outlet  temperatures  were  predicted  within  1.7  degrees  F  of 
that  reported.  The  outlet  temperatures  for  Run  20  were  predicted 
fairly  well  using  the  McAdams  correlation,  the  difference  in  shell 
outlet  temperature  as  given  by  experiment  and  theory  being  1.5  degrees  F, 
and  that  in  tube  outlet  temperature  being  1.8  degrees  F.  The  Stein  and 
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Begell  correlation  gave  temperature  predictions  within  the  region  of 
2.0  to  2.5  degrees  F.  The  poorest  agreement  resulted  for  Run  23,  in 
which  the  theoretical  heat  transfer  was  approximately  75  per  cent  of 
that  given  by  experiment.  It  is  thought  significant  that  this  run  was 
conducted  with  the  lowest  Reynolds  numbers  in  the  tube  and  shell,  both 
being  approximately  5000.  It  will  be  recalled  that  the  heat  transfer 
correlations  used  in  the  model  are  recommended  for  flow  conditions 
for  which  the  Reynolds  number  is  greater  than  10,000. 


5.3  Simulation  Study  using  the  Experimental  Data  obtained  from  the 

Test  Heat  Exchanger 

On  the  basis  of  the  previous  work  it  was  decided  to  use  McAdams' 
correlation  for  the  tube-to-annulus  film  coefficient  for  the  remainder 
of  the  simulation  studies.  The  test  heat  exchanger,  employed  in  later 
work  for  on-line  pulse  testing,  was  used  to  further  test  the  CSMP  model. 
The  basis  of  the  testing  was  not  only  the  model's  ability  to  predict 
steady  state  conditions,  as  in  the  previous  section,  but  dynamic  condi¬ 
tions  as  well.  To  obtain  the  necessary  data,  advantage  was  taken  of 
the  fact  that  the  sensing  and  control  devices  were  interfaced  with  the 
IBM  1800  computer.  A  series  of  process  programs,  explained  fully  in 
Chapter  15,  were  written  to  attain  three  primary  functions  -  data 
acquisition,  data  processing,  and  convenient  data  presentation. 

In  the  first  program,  WLHBL,  the  two  flows  and  four  temperatures 
were  recorded  and  logged  every  second,  for  a  user-specified  interval  of 
time.  These  data  were  then  averaged  and  a  heat  balance  calculation 
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performed;  the  average  values  and  the  result  of  the  heat  balance  were 
printed  out  on  one  of  the  DACS  Centre’s  typewriters.  The  process 
program  WLHBL  was  written  such  that,  after  initialization,  the  data 
logging,  averaging  and  processing  were  repeated  continuously  until  the 
user  decided  either  to  abort  or  to  continue  to  the  next  process  program. 
Either  decision  was  effected  by  appropriate  use  of  the  computer  console 
data  switches.  If  only  steady  state  runs  were  required,  the  user  would 
abort  at  this  stage  of  the  run. 

By  linking  to  the  next  and  subsequent  programs,  the  user  could 
obtain  dynamic  data.  Via  a  keyboard  request  feature,  the  user  is  able 
to  specify  a  pulse  shape  and  width  with  which  he  desires  to  force  the 
tube  flow  control  valve.  On  implementation  of  the  disturbance,  the 
computer  is  programmed  to  automatically  read  and  store  the  tube  flow 
and  outlet  shell  temperature  response  for  a  period  of  forty  seconds. 

By  limiting  the  pulse  duration  to  under  ten  seconds,  enough  time  is 
allowed  for  the  transients  to  die  out.  The  integer  data  collected 
during  the  run  is  converted  to  engineering  values  (U.S.  gallons  per 
minute,  and  degrees  Fahrenheit)  by  appropriate  conversion  programs. 

The  processed  data  may  then  be  plotted  out  on  the  oscilloscope  or  X-Y 
plotter,  or  it  may  be  punched  out  on  computer  cards.  One,  two  or  all 
three  of  these  options  may  be  specified  by  the  user  via  the  on-line 
keyboard. 

For  the  steady  state  simulation,  experimental  data  for  nine 
different  runs  were  obtained;  the  tube  flows  were  6,  9  and  12  U.S. 
gallons  per  minute,  the  shell  flows  were  15,  20  and  25  U.S.  gallons  per 
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minute,  and  the  inlet  fluid  temperatures  were  in  the  region  of 

64.0  degrees  F  and  147.0  degrees  F  for  the  tube  and  shell  respectively. 

The  results  of  the  steady  state  simulation  runs  appear  in 
Table  3.  As  shown,  the  ranges  of  tube  liquid  and  annulus  liquid 
Reynolds  numbers  were  26,000  to  50,000,  and  42,000  to  71,000 
respectively.  It  is  believed  that  these  flows  are  sufficiently  large  in 
range  to  thoroughly  test  the  model's  steady  state  temperature  predictions. 

All  of  the  nine  runs  predicted  a  heat  transfer  rate  of  within 
+  7.0  per  cent  of  that  actually  recorded.  Because  of  the  heat  loss 
from  the  shell,  the  experimental  results  all  indicated  a  Q  /Q  ratio 

u  -L 

greater  than  unity.  If  this  heat  loss  is  taken  into  account,  comparison 
of  the  actual  and  predicted  heat  transfer  rates  range  between  -  1.5  per 
cent  to  +  7.0  per  cent.  The  majority  of  the  runs  still  indicate, 
however,  that  the  model  predicts  lower  heat  transfer  rates  than  those 
experienced  in  practice.  There  are  baffles  in  the  test  heat  exchanger 
which  are  used  to  ensure  correct  spacing  between  the  tube  and  shell 
along  its  10-foot  length.  These  baffles  result  in  higher  turbulence 
levels  and  higher  values  of  the  heat  transfer  film  coefficients  than 
those  calculated  for  the  ideal  case.  This  is  suspected  to  be  the  prime 
reason  for  the  discrepancy  in  the  heat  transfer  rates.  Another 
possibility  is  that  the  heat  transfer  correlations,  for  this  range  of 
Reynolds  numbers,  do  in  fact  give  lower  than  actual  values  for  the  film 
coefficients . 

Simulation  of  the  experimental  heat  exchanger  was  more  success¬ 
ful  than  the  simulation  of  Fisher's  data.  It  was  believed  that  the 
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main  reason  for  this  was  that  the  Reynolds  numbers  attained  in  the 
test  heat  exchanger  studies  were  greater  than  10,000,  whereas  for  the 
particular  runs  chosen  from  Fisher’s  data,  the  Reynolds  numbers  were 
between  5,000  and  7,000. 

Considering  the  limited  accuracy  of  the  thermocouples  (+0.5 
degrees  F) ,  and  the  heat  transfer  correlations  (+  10  per  cent),  the 
results  of  the  steady  state  simulation  work  are  very  satisfactory. 

For  the  dynamic  simulation  studies,  three  experimental  runs 
were  made  in  which  the  forcing  function  (tube  flow),  and  response 
function  (outlet  shell  temperature)  were  logged  and  plotted.  Two  of 
the  runs  employed  a  step;  the  remaining,  a  rectangular  pulse,  to  excite 
the  system.  The  pertinent  conditions  for  each  run  are  given  in  Table  4; 
the  forcing  functions  and  the  resulting  transients  are  given  in  Table  5 
and  shown  in  Figures  3  to  8 .  The  experimental  runs  were  simulated,  and 
the  temperature  responses  as  given  by  the  CSMP  model  are  shown  as  the 
dotted  lines  in  Figures  4,  6  and  8.  For  all  three  runs,  the  actual  and 
theoretical  responses  are  dynamically  similar.  However,  the  experi¬ 
mental  results,  unlike  those  of  the  model,  indicate  a  time  lag  of 
approximately  one  second  before  the  transient  occurs.  This  is  partly 
the  result  of  the  dynamic  lags  of  the  elements  in  the  control  and 
sensing  loops.  Two  main  contributors  to  this  time  delay  are  the  control 
valve  and  thermopile. 

Ignoring  the  time  delay,  the  transient  responses,  actual  as  well 
as  theoretical,  appear  to  exhibit  those  of  a  simple  first-order  system; 
this  concurs  with  the  results  of  previous  investigators  of  heat 
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exchanger  dynamics  (56,  26,  42).  On  the  basis  that,  for  a  first-order 
system,  63.2  per  cent  of  the  final  change  in  the  response  variable  to 
a  step  forcing  function  occurs  after  one  time  constant,  Runs  10  and  11 
were  analyzed.  The  time  constant  given,  both  by  experiment  and  theory, 
is  approximately  2.5  seconds. 
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TABLE  5 


Shell  Outlet  Temperature  Responses 
as  given  by  the  CSMP  Model 


Run 

#10 

Run  #11 

Run  #12 

Time 

(secs) 

Tube  Flow 
(ft/ sec) 

Shell  Temp. 
(°F) 

Tube  Flow 
(ft/ sec) 

Shell  Temp . 
(°F) 

Tube  Flow 
(ft/ sec) 

Shell  Temp. 
(°F) 

0 

6.15 

130.75 

3.78 

132.27 

4.92 

148.70 

1 

6.15 

130.75 

3.78 

132.27 

4.92 

148.70 

2 

6.15 

130.75 

3.78 

132.27 

4.92 

148.70 

3 

6.15 

130.75 

3.78 

132.27 

4.92 

148.70 

4 

6.15 

130.75 

3.78 

132.27 

4.92 

148.70 

5 

6.15 

130.75 

3.78 

132.27 

4.92 

148.70 

6 

4.00 

131.50 

6.33 

131.42 

3.20 

149.60 

7 

4.00 

132.17 

6.33 

130.61 

3.20 

150.40 

8 

4.00 

132.52 

6.33 

130.14 

3.20 

150.82 

9 

4.00 

132.70 

6.33 

129.88 

3.20 

151.02 

10 

4.00 

132.78 

6.33 

129.76 

3.20 

151.12 

11 

4.00 

132.81 

6.33 

129.71 

4.92 

150.26 

12 

4.00 

132.83 

6.33 

129.68 

4.92 

149.47 

13 

4.00 

132.84 

6.33 

129.67 

4.92 

149.06 

14 

4.00 

132.84 

6.33 

129.67 

4.92 

148.86 

15 

4.00 

132.84 

6.33 

129.66 

4.92 

148.77 

16 

4.00 

132.84 

6.33 

129.66 

4.92 

148.73 

17 

4.00 

132.84 

6.33 

129.66 

4.92 

148.71 

18 

4.00 

132.84 

6.33 

129.66 

4.92 

148.71 

19 

4.00 

132.84 

6.33 

129.66 

4.92 

148.70 

20 

4.00 

132.84 

6.33 

129.66 

4.92 

148.70 
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FIG.  3  TUBE  FLOW  FORCING  FUNCTION  -  RUN  10 


EXPERIMENTAL 
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FIG.  4  COMPARISON  OF  EXPERIMENTAL  AND  THEORETICALLY  PREDICTED 
SHELL  FLUID  OUTLET  TEMPERATURE  TRANSIENTS  -  RUN  10 
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FIG.  5  TUBE  FLOW  FORCING  FUNCTION  -  RUN  11 
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FIG.  6  COMPARISON  OF  EXPERIMENTAL  AND  THEORETICALLY  PREDICTED 
SHELL  FLUID  OUTLET  TEMPERATURE  TRANSIENTS  -  RUN  11 
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FIG.  7  TUBE  FLOW  FORCING  FUNCTION  -  RUN  12 
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FIG.  8  COMPARISON  OF  EXPERIMENTAL  AND  THEORETICALLY  PREDICTED 
SHELL  FLUID  OUTLET  TEMPERATURE  TRANSIENTS  -  RUN  12 


6  CONCLUSIONS 


The  steady  state  and  dynamic  simulation  results  are  within  the 
limits  of  accuracy  imposed  by  the  sensing  devices  and  heat 
transfer  correlations. 

The  success  of  the  simulation  studies  was  thought  to  be  due  to 
two  prime  factors.  First,  a  complete  correlation  for  each  of 
the  three  heat  transfer  coefficients  was  incorporated  in  the 
model.  Secondly,  the  effect  of  temperature  on  the  thermal 
conductivity  and  viscosity  of  water  was  not  ignored  when  calculat¬ 
ing  the  film  coefficients. 

The  heat  exchanger  appears  to  exhibit  the  dynamics  of  a 
transportation  lag  and  first-order  system  in  series.  The 
transportation  lag  was  about  1.0  second,  and  the  value  of  the 
first-order  system  time  constant  approximately  2.5  seconds. 

Privott  (46),  in  his  investigation  of  heat  exchanger  dynamics, 
posed  the  question  as  to  whether  heat  transfer  film  coefficient 
correlations,  the  derivations  of  which  are  based  on  steady  state 
studies,  are  valid  for  unsteady  state  conditions.  The  ability 
of  the  model  used  in  this  study  to  predict  the  shell  fluid 
outlet  temperature  transients  indicates  that  the  correlations 
are  indeed  satisfactory  for  the  dynamic  case. 

A  satisfactory  simulation  of  the  10-foot  long  heat  exchanger  was 
realized  using  only  a  four-lump  model.  It  is  not  known  what 
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value  of  the  heat  exchanger  length- to-lump  ratio  has  to  be 
reached  before  a  detrimental  effect  occurs  in  the  simulation. 


A  digital  simulation  language  approach,  although  not  simple, 
provides  a  convenient  means  of  solving  distributed  parameter 
problems.  As  to  the  worth  of  this  method,  it  should  be 
recalled  that  prior  to  this  study  all  previous  workers  on  heat 
exchanger  dynamics  have  resorted  to  simplifying  expressions 
for  the  film  coefficients  before  solving  their  models.  Also,  to 
solve  the  general  four-equation  model  on  an  analog  computer 
would  impose  a  formidable,  if  not  impossible,  task. 
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7  INTRODUCTION  TO  PULSE  TESTING 

7.1  Introduction 

The  engineer  concerned  with  designing  a  control  system  for  a 
chemical  process  or  plant  has  at  hand  a  wealth  of  theory  on  sophisticated 
design  techniques.  The  majority  of  these  techniques  have  been  evolved 
on  two  assumptions.  Firstly,  the  process  is  a  linear  system,  which 
means  that  the  relation  between  an  output  response  from  an  input 
disturbance  may  be  described  with  a  linear  integro-dif f erential  equation 
with  constant  coefficients.  The  second  assumption  is  that  the 
mathematical  model  or  transfer  function  for  the  system  is  known.  For 
the  electrical  engineer  who  deals  with  components,  the  process 
characteristics  of  which  are  usually  well  defined,  this  latter  assump¬ 
tion  may  be  valid.  However,  for  the  chemical  engineer  this  is  rarely 
the  case;  the  complexities  of  most  chemical  unit  operations  tend  to 
defy  simple  a  priori  analysis.  This  factor  has  tended  to  deter  engineers 
in  the  chemical  industry  from  using  the  process  dynamic  approach  to 
design,  and  to  rely  simply  on  steady  state  information. 

Hougen  (28)  states  that  more  interest  is  being  shown  in  process 
characterization  and  control  than  ever  before,  but  even  so,  "surprisingly 
few  in  chemical  engineering  are  engaged  in  procuring  and  utilizing 
dynamic  information".  The  obvious  benefit  of  having  reliable  dynamic 
information  is  that  the  sophisticated  concepts  of  process  control  may 
be  applied.  Perhaps  not  just  simple  feedback,  but  feed-forward  and 
adaptive  control  schemes  may  be  implemented.  Also,  knowing  more  about 
the  process  results  in  better  design  and,  possibly,  better  plant 
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scale-up.  Other  benefits  accruing  from  knowledge  of  plant  dynamics 
and  perhaps  not  so  obvious  are  that  it  could  be  useful  in  the  deter¬ 
mination  of  kinetic  mechanisms,  and  it  also  could  help  in  diagnosing 
malfunctions  and  improper  design  of  chemical  equipment. 


7.2  Methods  of  Obtaining  a  Model 

Five  general  methods  for  obtaining  a  model  of  a  process  are: 

(i)  Analytic 

(ii)  Frequency  response 

(iii)  Step  forcing 

(iv)  Random  inputs 

(v)  Pulse  inputs. 

The  first  is  the  theoretical  approach,  the  latter  four  are  experimental 
techniques  which  are  classified  by  the  manner  in  which  the  system 
under  study  is  excited. 

(i)  Analytic  Method: 

The  nature  of  most  chemical  processes  does  not  lend  it¬ 
self  to  simple  dynamic  analysis.  Heat,  mass,  momentum 
and  chemical  reaction  may  be  inherent  in  a  process,  and 
a  perusal  of  the  book  written  by  Franks  (24)  will  indi¬ 
cate  to  the  reader  the  large  number  of  equations  which 
are  needed  to  define  even  the  simplest  of  systems.  The 
result  is  that,  for  most  chemical  processes,  the  engineer 
tends  to  rely  only  on  steady  state  information  for  design 
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(ii)  Frequency  Response  Method: 

The  classical  experimental  approach  of  obtaining  dynamic  infor¬ 
mation  is  by  frequency  response  testing.  An  input  variable 
is  varied  sinusoidally,  time  is  allowed  for  transients  to  die 
out,  and  the  sinusoidally  varying  output  variable  is  recorded. 
The  amplitude  ratio  and  phase  lag  between  the  two  sinusoids 
at  that  frequency  is  calculated.  This  procedure  is  repeated 
for  a  whole  range  of  frequencies  and  the  information  con¬ 
veniently  expressed  in  the  form  of  a  Bode  diagram.  Because 
the  method  readily  lends  itself  to  filtering  and  signal 
averaging  techniques,  and  because  it  is  easy  to  check  non- 
linearities  in  the  process  by  simple  observation  of  the  input 
and  output  sinusoids,  high  confidence  may  be  expected  in  the 
accuracy  of  the  resulting  frequency  information.  There  are, 
however,  drawbacks  when  trying  to  implement  frequency  response 
testing  in  a  plant.  The  method  is  exceedingly  time  consuming, 
it  may  be  difficult  to  realize  a  good  input  sinusoid,  and 
the  pecuniary  loss  due  to  off-specification  product  may  be 
severe . 

(iii)  Step  Response  Method: 

As  the  name  indicates,  a  step  is  introduced  into  the  system; 
the  time  histories  of  the  forcing  and  response  variable  are 
recorded.  By  appropriate  mathematical  manipulation,  the 
frequency  response  of  the  system  may  be  obtained.  There 
are  a  number  of  disadvantages  to  this  method,  the  most 
serious  being  the  loss  of  accuracy  at  high  frequencies. 
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(iv)  Random  Input  Method: 

By  using  auto-  and  cross-correlation  methods  it  is 
theoretically  possible  to  derive  the  frequency  response 
of  a  system  by  using  pure  random  signals.  It  has  been 
applied  successfully  (30),  although  a  serious  limitation 
of  this  method  is  that  a  large  number  of  data  points  are 
required  in  order  to  obtain  any  useful  information;  also, 
sensing  devices  have  to  be  very  accurate. 

(v)  Pulse  Testing  Method: 

Any  waveform  may  be  represented  as  the  summation  of  a 
series  of  sinusoids  of  different  amplitude.  Thus,  if  an 
input  variable  to  a  system  is  pulsed,  the  dynamics  of  the 
process  are  excited  by  all  the  component  frequencies  of 
that  pulse.  This  means  that,  from  a  well-planned  pulse 
test,  it  is  possible  to  extract  the  same  information  as 
a  whole  series  of  frequency  response  tests.  In  Part  VI 
of  their  Chemical  Engineering  Refresher  Series,  Murrill, 
Pike  and  Smith  (44)  state,  "Of  all  the  methods  presented 
in  this  series  for  determining  dynamic  models  of 
machinery,  process  equipment,  reactors,  etc.,  pulse 
testing  techniques  are  probably  the  ones  with  the- 
greatest  potential". 

In  this  part  of  the  thesis  theoretical  consideration  is  given  to 
various  aspects  of  pulse  testing.  The  study  culminates  in  Chapter  13 
in  which  the  reader  is  given  a  set  of  guidelines  for  designing  and 
analyzing  the  results  of  a  pulse  testing  experiment. 
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8  PULSE  TESTING  LITERATURE  REVIEW 


Pulse  testing  seems  to  have  had  its  origin  in  the  aircraft 
industry.  In  the  late  1940’s  the  appearance  of  transonic  and  super¬ 
sonic  aircraft  increased  the  problems  of  designing  satisfactory 
automatic  and  stabilizing  control  systems.  Theoretical  analsyis 
alone  of  the  performance  of  an  aircraft  was  no  longer  felt  adequate 
and  it  was  necessary,  therefore,  to  resort  to  flight-test  procedures 
to  substantiate  the  design  theory. 

Seamens,  Bromberg  and  Payne  (49)  first  described  the  concept  of 
a  performance  function  relating  the  output  response  of  a  system  to  an 
input  disturbance.  Seamens,  Blasingame  and  Clementson  (48),  on  the 
advice  of  the  latter  author,  actually  applied  pulse  forcing  instead  of 
the  usual  step  forcing  to  obtain  certain  transfer  functions  of  a 
U.S.  Air  Force  B52  bomber.  These  authors  described  incremental  step 
and  triangular  area  approximations  for  solving  the  integrals  arising 
from  the  resulting  data  analysis. 

Walters  and  Rea  (57)  also  demonstrated,  about  this  time,  the 
possibility  of  determining  the  frequency  response  characteristics  of 
a  system  from  an  arbitrary  input  signal.  They  described  a  matrix 
multiplication  method  for  evaluating  the  Fourier  coefficient  integrals. 
Bollay  (4)  gives  a  general  description  of  the  progress  up  to  1951  of 
determining  aircraft  dynamics  using  various  methods;  included  is  a 
review  of  pulse  testing.  Smith  and  Triplett-  (52)  described  methods 
and  techniques  for  obtaining  airplane  frequency  response  characteris¬ 
tics  from  transient  flight  measurements.  The  authors  discussed  the 


. 
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effect  of  frequency  content  of  the  input  pulse  and  the  relationships 
of  pulse  shape  and  spectral  content. 

In  1954,  Draper,  McKay  and  Lees  (16)  provided  the  earliest 
authoritative  reference  on  pulse  testing.  The  mathematical  background 
necessary  for  the  numerical  evaluation  of  Fourier  transforms  is  dis¬ 
cussed  in  detail.  In  addition,  the  results  of  applying  two  different 
pulse  shapes,  rectangular  and  displaced  cosine,  to  a  first-  and  second- 
order  system  are  given.  Some  of  the  same  material  appears  in  a  paper 
by  Lees  (36).  This  author  presents  a  broad  view  of  the  field  of 
process  identification  including  background  concepts  necessary  for  the 
interpretation  of  dynamic  data,  as  well  as  criteria  for  selecting  a 
test  method. 

Lees  and  Hougen  (38)  in  1956  published  a  paper  in  which  they 
compared  the  pulse  and  direct  frequency  methods  applied  to  a  shell  and 
tube  heat  exchanger.  They  concluded  that  the  cost  of  pulse  testing 
should  be  considerably  less  than  direct  sinusoidal  forcing  methods. 

This  work  appears  to  be  the  first  to  herald  the  introduction  of  pulse 
testing  to  the  chemical  industry.  Their  work  was  later  extended  by 
Morris  (42).  Hougen  and  Walsh  (29)  described  other  applications  of 
pulse  testing  pertinent  to  chemical  engineering.  Head,  Hougen  and 
Walsh  (25)  used  pulse  testing  information  to  verify  a  mathematical 
model  describing  the  dispersion  of  a  solute  in  a  flowing  fluid  system. 
Another  paper,  also  concerned  with  a  mass  transfer  operation,  is  by 
Marino,  Perna  and  Stutzman  (39).  These  authors  applied  both  pulse 
forcing  and  direct  sinusoidal  forcing  of  the  liquid  reflux  return  of 
a  24-plate  distillation  column,  and  measured  the  resulting  temperature 
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transients  at  various  plate  locations.  From  the  resulting  frequency 
diagrams  it  was  concluded  that  the  system  could  be  approximated  by 
linear  first-order  equations. 

Following  the  work  of  Smith  and  Triplett  (52) ,  a  number  of 
workers  have  been  concerned  with  the  actual  study  of  pulse  testing  as 
a  process  identification  technique  per  se.  Dreifke  (17)  made  an 
extensive  study  of  pulse  shapes  and  widths,  and  the  resulting  accuracy 
of  the  dynamic  analysis.  Some  of  the  conclusions  of  Dreifke' s  work 
appear  in  a  Monograph  written  by  Hougen  (28).  Although  this  Monograph 
deals  with  the  subject  of  process  dynamics  in  general,  a  major  portion 
of  it  is  devoted  to  pulse  testing.  Clements  (8)  investigated  the 
theoretical  as  well  as  practical  problems  of  the  subject,  and  presented 
some  interesting  conclusions.  His  experimental  work  involved  the 
application  of  pulse  testing  to  the  study  of  the  dynamics  of  an 
extraction  column.  A  report  by  Lees  and  Dougherty  (37)  discusses  the 
errors  involved  in  the  numerical  computation  of  the  Fourier  transform 
of  the  input  and  response  functions,  and  one  of  their  conclusions  is 
that  the  upper  limit  for  valid  frequency  data  is  o>  At  £  O.Itt,  not  u 
as  predicted  by  the  sampled  data  theorem. 

In  addition  to  the  above,  reviews  dealing  with  the  theory  and 
application  of  pulse  testing  have  been  presented  by  Clements  and 
Schnelle  (9),  and  Dreifke  and  Hougen  (18). 
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9  MATHEMATICAL  BACKGROUND 


The  mathematical  background  of  pulse  testing  is  well  documented 
in  the  literature;  the  approach  in  this  thesis  follows  that  of 
Clements  (8).  In  essence,  the  following  mathematical  outline  details 
how  pulse  testing  data,  obtained  in  the  time  domain,  may  be  processed 
to  give  the  input-output  relationship  of  a  system  in  the  frequency 


domain. 


The  transfer  function,  G(s) ,  of  a  system,  is  given  by  the  ratio 
of  the  Laplace  transforms  of  its  output  and  input  time  functions: 


. ..  9-1 


The  upper  limits  of  the  integrals  are  infinity,  but  in  actual 
practice  the  pulse  duration  and  the  transient  response  will  have  finite 


time  durations  given  by  0  <  t <  T  ,  and  0  <  t  <  T  respectively. 

x  y 


Equation  9-1  may  be  further  modified  by  replacing  the  complex  variable 
of  the  Laplace  transform,  s,  by  ju),  resulting  in  the  transfer  function 
being  expressed  in  the  frequency  domain  (5): 


. ..  9-2 


The  equation  is  simplified  by  using  the  identity, 


e 


-jut 


cos  wt  -  j  sin  03t 


G(jco) 


_  f Q  y(t)  cos  cot  dt  -  j  f  ■'y(t)  sin  cot  dt 

~T  T 

^  ■** 

fo  x(t)  cos  cot  dt  -  j  /Qxx(t)  sin  cot  dt  ...  9-3 


Defining  the  following  relationships, 

T 

A  =  /  yy(t)  cos  cot  dt 

o  J 

T 

B  =  S J y(t)  sin  cot  dt 
T 

C  =  /  Xx(t)  cos  cot  dt 

0 

T 

D  =  /  Xx(t)  sin  cot  dt 

o 


. ..  9-4 


. . .  9-5 


. ..  9-6 


. ..  9-7 


the  frequency  response  may  be  conveniently  expressed  as, 


G(jco)  =  A  -  jB 
C  -  jD 


.  ..  9-8 


Multiplying  the  top  and  bottom  of  the  equation  by  the  complex  conjugate 
of  the  denominator,  one  can  obtain  the  real  and  imaginary  parts  of  G(jco) 


Re  (co)  =  AC  +  BD 
C2  +  D2 


. ..  9-9 


Im  (co)  =  AD  -  BC 
C2  +  D2 


.  ..  9-10 


It  thus  follows  that  the  magnitude  ratio  and  phase  angle  for  any  given 
frequency  is  given  by. 


■  '•  ■  V 


■ 
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M.R.  (to)  =  /[Re(oo)  ] 2  +  [Im(m)]2 


.  ..  9-11 


and 


<f>(co) 


=  tan  Im(m) 
Re(w) 


. ..  9-12 


In  pulse  testing  work,  it  is  important  to  know  the  magnitude  of  the 
Fourier  amplitudes  of  the  forcing  function.  The  frequency  content,  as 
it  is  termed,  of  a  pulse,  is  given  by  the  absolute  magnitude  of  its 
Fourier  transform, 


s(o)> 


/xx(t)  e-Jut  dt 


=  ^C2  +  D2 


. ..  9-13 


It  is  convenient,  when  discussing  the  frequency  content  of  different 
pulses,  to  normalize  s(to)  by  dividing  it  by  the  content  at  zero  frequency, 


sMn 


sM 

s(0) 


_ s(co) _ 

[total  area  under  x(t)] 


s(aQ 

T 

/  xx(t)  dt 

o 


9-14 
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Using  the  above  relationships,  pulse  test  time  data  may  be 
reduced  to  frequency  response  data.  The  major  numerical  problem  is 
the  evaluation  of  the  real  and  imaginary  parts  of  the  Fourier 
transforms,  i.e.,  the  quadrature  products  A,  B,  C  and  D;  this  problem 
is  discussed  in  detail  in  Chapter  11. 
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10  VALIDITY  OF  PULSE  TESTING  RESULTS 


10.1  Theoretical  Study  of  the  Causes  of  Inaccuracies  occurring  in 

Pulse  Testing  Results 

Experimental  and  theoretical  frequency  response  data,  determined 
from  pulse  testing  results,  may  suffer  from  inaccuracies  due  to  one  or 
more  of  the  following: 

(i)  Nonlinearities  in  the  process 

(ii)  Frequency  content  of  the  forcing  function  as  it  approaches 
zero 

(iii)  Aliasing  of  frequencies 

(iv)  Overlapping  of  frequencies  due  to  a  non-bandwidth  limited 
forcing  function 

(v)  Numerical  approximation  errors 

(vi)  Truncation  of  the  time  data. 

(i)  Nonlinear  process: 

If  a  pulse  is  applied  to  a  linear  process,  then  the 
sinusoidal  components  of  the  pulse  will  result  in  sinu¬ 
soids  emanating  from  the  system  with  identical  frequen¬ 
cies,  although  differing  in  magnitude  and  phase.  If  the 
system  is  nonlinear,  the  output  wave  forms  will  be 
distorted,  giving  rise  to  corresponding  distortions  in 
the  frequency  diagram. 


r 


(ii)  Frequency  Content  of  the  Pulse: 


A  pulse  has  its  own  peculiar  frequency  content,  given  by 
the  following  expression: 


FT  [f (t) ] 


/  f(t)  e 

o 


-joot 


dt 


As  indicated  previously,  if  this  is  divided  by  the  area 

under  the  pulse,  the  normalized  frequency  content,  s(u>)  , 

n 

is  obtained.  The  normalized  frequency  content  curves  of 
pulses  used  by  the  author  in  the  experimental  work  of 
this  thesis  are  shown  in  Chapter  16. 

Note  that  the  normalized  frequency  content  approaches 
zero  as  oo  tends  to  infinity: 

s  ( co )  -*■  0  as  oo  -*  00 

n 

The  result  is  that  the  amplitudes  of  the  high  frequency 
components  of  the  pulse  are  too  small  to  excite  the 
system  sufficiently  to  extract  any  useful  information. 
Also,  as  the  quadrature  products  A,  B,  C  and  D  tend  to 
zero,  numerical  roundoff  errors  become  significant. 

Both  factors  tend  to  invalidate  the  results  at  higher 
frequencies . 

(iii)  Aliasing: 

Aliasing  results  from  treating  a  continuous  function  as 
a  sampled  data  series.  The  numerical  techniques  employed 
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to  reduce  pulse  testing  information  to  frequency 
response  results  require  that  the  input  data  be  expressed 
as  a  series  of  discrete  values  which  represent  the  con¬ 
tinuous  time  function.  The  data  may  be  logged  by  a 
digital  computer  at  equidistant  intervals  of  time,  or  the 
data  recorded  continuously  by  means  of  some  plotting 
device  and  the  discrete  values  read  from  this  record. 

Either  method  reduces  the  continuous  time  function  to  a 
discrete  set  of  values.  Suppose  the  sampling  interval 
is  At  seconds,  then  the  frequency  of  the  sine  wave  to 
pass  exactly  through  the  sampling  interval  points  is 
given  by  Tr/At  radians  per  second.  This  particular 
frequency  is  known  as  the  Nyquist  frequency.  It  is 
evident,  however,  that  sine  waves  of  higher  frequencies 
given  by 


oo,  =  kn_,  k  =  1,  2, 
k  At 


...  10-1 


will  also  pass  through  the  sampling  points,  and  as  far  as 
the  data  is  concerned,  the  Nyquist  frequency  and  its 
associated  "folding"  frequencies  are  indistinguishable. 
Since  there  is  no  information  about  x(t)  between  the 
sampling  points,  there  is  no  means  of  discerning  the 
amplitude  of  frequencies  higher  than  the  Nyquist 
frequency.  Thus,  sampled  data  analysis  would  therefore 
warn  investigators  of  assuming  their  frequency  data  is 


valid  past  the  Nyquist  frequency, 


oo-  .  =oo  =  it  rads/sec. 

llmlt  n  At 


...  10-2 


•  - 


' 
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(iv)  Overlapping  of  Frequencies: 

The  sampling  process  inherent  in  digital  data  acquisition 
systems  may  be  represented  by  a  pulse-amplitude  modula¬ 
tion  process  that  converts  a  continuous  signal  into  a 
train  of  amplitude  pulses.  For  a  zero  order  hold  sampler, 
the  continuous  function  f(t)  is  represented  by  the 
stepped  curve  f*(t),  as  shown  in  Figure  9  below: 


Fig.  9  Stepped  Curve  Representation  of  a  Continuous 

Function  as  given  by  a  Zero  Order  Hold  Sampler 


By  inspection  of  the  two  responses  one  may  expect  to 
find  much  higher  frequency  components  in  f*(t)  than  in 
f(t)  because  of  the  discontinuity  characteristic  of  the 
pulse  train.  One  can  describe  the  operation  of  the  pulse- 
amplitude  modulator  by 

f*(t)  =  p(t)  f(t)  ...  10-3 


where  f(t)  is  the  continuous  input  signal,  f*(t)  is  the 


. 
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output  signal  of  the  sampler,  and  p(t)  denotes  the  unit- 
pulse  train  carrier.  If  the  sampling  period  T  is  constant, 
p(t)  is  a  per iodicf unction  and  can  be  expanded  into  a 
Fourier  series: 


p(t) 


. ..  10-4 


where  ms  is  the  sampling  frequency  in  radians  per  second 
and  u)s  is  related  to  the  sampling  period  T  by 


oos  =  2tt/T 


. ..  10-5 


The  output  of  the  sampler  may  now  be  written  as 


00 

f*(t)  =  z  C  f(t)  e^na)st  ...  10-6 

n  =  -  00  n 


By  using  the  shifting  theorem  of  the  Fourier  transform 
which  states  that 


cfjf(t)  e^naJst]  =  F(ja>  -  jna>s)  ...  10-7 


it  follows  that 


C  F 
n 


(jw  -  jnws) 


F*(ju>) 


•  •  • 


10-8 


(  it  “  =  *'»  ( : 
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The  latter  equation  may  also  be  written  as 


00 

F*( ja>)  =  n=z_oo  Cn  F  (ju>  +  jnoos)  ...  io-9 


since  n  is  summed  from  to  +00. 
The  magnitude  of  F*(jaj)  is  simply 


F*(ju) 

— 

oo 

l  C  F  (jo)  +  jnoos) 

n  =  -  co  n 

Thus  F*(jw)  consists  not  only  of  the  fundamental  compon¬ 
ent  of  F(jw),  but  also  an  infinite  number  of  spectra 
called  the  complimentary  components. 

Consider  a  band  limited  function  which  has  a  frequency 
content  of  zero  at  00  ,  also  that  the  sampling  of  this 

r 

function  occurs  at  more  than  twice  co  .  The  amplitude 
spectra  as  given  by  Equation  10-10  will  be  as  shown  in 
Figure  10.  The  original  signal  can  be  recovered  from 
the  spectrum  by  an  ideal  low  pass  filter.  Suppose,  how¬ 
ever,  that  co  <  2oo  ;  as  Figure  11  indicates,  the  spurious 

S  i* 

or  complimentary  spectra  overlap  the  fundamental  compon¬ 
ent,  and  distortion  of  F*(jw)  will  occur,  and  the  input 
signal  cannot  be  reconstructed  even  by  means  of  an 


ideal  filter. 
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Fig.  10  Amplitude  Spectra  of  Band  Limited  Function 


Fig.  11  Spectra  Distortion  caused  by  Overlapping  of 
Frequencies 


hbIdi  ><  -  f' 
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In  reality,  there  is  no  such  thing  as  a  band  limited 
signal,  thus,  even  with  application  of  a  low  pass  filter, 
and  a  sampling  rate  ws  >  2toF,  distortion  of  F*(jw)  will 
still  occur. 

(v)  Numerical  Approximation  Errors: 

Lees  and  Dougherty  (37)  investigated  the  inaccuracies 
of  numerically  calculating  the  Fourier  transform  of  the 
input  and  response  resulting  from  a  pulse  test.  They 
were  particularly  concerned  with  how  well  the  summation 
factor 


Zf(i  At)  exp  (-jcoiAt)  where  i  =  1,  . .  N 

represents  the  ideal  transform 

/f(t)  e"jaJt  dt. 

They  concluded  that  the  approximation  is  reasonably  good 
up  to  u)At  =  O.Itt.  This  indicates  that  at  least  ten  times 
as  many  data  points  are  required  per  interval  as  predicted 
by  the  sample  data  theorem. 

(vi)  Truncation  of  Time  Data: 

This  problem  was  studied  by  Dreifke  (17),  who  concluded 
that  "moderate  truncation  of  the  time  response,  say  up 
to  23  per  cent  or  50  per  cent  of  the  time  history,  has 


f‘ 


negligible  deleterious  effects  on  the  related  experi¬ 
mental  frequency  response  curve". 
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10.2  Summary  and  Conclusions  of  Study 

The  usual  purpose  of  conducting  a  pulse  test  is  to  provide 
frequency  response  information  from  which  the  engineer  might  fit  a 
model  to  the  process.  Any  distortions  of  the  amplitude  ratio  and 
phase  angle  curves  due  to  nonlinearities  in  the  process  tend  to  be 
"smoothed  out"  in  the  curve  fitting  procedure,  and  therefore  are 
usually  of  no  serious  concern.  Because  it  is  not  feasible  to  have  a 
bandwidth  limited  signal  in  reality,  distortions  arise  due  to  over¬ 
lapping  frequencies.  For  the  same  reason  given  above,  this  usually 
presents  no  serious  problem  to  the  chemical  engineer.  This  is  not  the 
case,  however,  in  the  field  of  electronics,  where  the  "unscrambling" 
of  frequency  data  by  spectral  smoothing  techniques  is  often  an 
essential  operation  (12). 

Rarely  should  the  engineer  encounter  a  situation  where  the  time 
data  has  to  be  truncated;  even  so,  from  the  findings  of  Dreifke,  a  severe 
cutoff  of  25  per  cent  to  50  per  cent  of  the  time  record  should  not 
seriously  affect  the  frequency  results. 

The  problem  of  aliasing  is  overshadowed  by  the  more  restrictive 
criteria  of  Lees  and  Dougherty  (37)  pertaining  to  numerical  approximation 
errors.  It  will  be  recalled  that  these  authors  recommend  that  ten  times 
the  number  of  time  data  points  be  used  as  given  by  the  sampled  data  theorem. 
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The  normalized  frequency  content  of  the  input  pulse  is  a  very 
important  guide  in  determining  the  maximum  frequency  to  which  the 
computed  frequency  response  curves  may  be  assumed  valid.  Clements  (8), 
on  the  basis  of  the  results  of  previous  workers  and  on  his  own  studies, 
stated  that  deviation  in  the  frequency  curve  usually  sets  in  when  the 
normalized  frequency  content  drops  to  a  value  of  approximately  0.3.  It 
was  also  stated  that  it  is  sometimes  possible  to  recover  reliable 
frequency  response  information  between  the  zeros  of  the  frequency 
content  curve  but,  since  this  is  the  exception  rather  than  the  rule  in 
most  practical  situations,  this  is  not  a  recommended  practice. 

In  summary,  for  the  chemical  engineer,  the  two  most  important 
factors  regarding  the  validity  of  pulse  testing  results  are  usually 
frequency  content  and  numerical  approximation  problems.  Fortunately, 
numerical  approximation  errors  may  be  avoided  at  the  design  stage  of 
a  pulse  test  by  simply  applying  the  Lees  and  Dougherty  criteria. 
Avoidance  of  errors  arising  from  unfavourable  frequency  content  curves 
is,  however,  not  such  an  easy  problem  to  solve.  In  essence,  for  a 
pulse  test  conducted  in  a  plant,  the  exact  pulse  shape,  and  therefore 
frequency  content  curve,  is  usually  not  known  prior  to  the  testing. 
Since  the  exact  frequency  content  of  the  pulse  can  only  be  calculated 
on  completion  of  the  plant  test  it  is  therefore  limited  as  a  design 
tool.  As  already  indicated,  however,  it  is  of  prime  importance  when 
assessing  the  validity  of  the  frequency  results. 


- 
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11  COMPARISON  OF  DATA  ANALYSIS  METHODS 

11.1  Introduction 

In  Chapter  9,  it  was  shown  that  the  transfer  function  G(jco) 
for  a  system  could  be  expressed  as  the  ratio  of  the  Fourier  transforms 
of  the  system  output  and  input.  Further,  this  expression  could  be 
broken  down  into  the  four  quadrature  products: 

T 

A  =  f  y  y(t)  cos  cot  dt 

o 

T 

B  =  f  ^  y(t)  sin  cot  dt 
T 

C  =  /  X  x(t)  cos  cot  dt 


T 

D  =  /  X  x(t)  sin  cot  dt 

o 

The  question  now  arises  as  to  which  numerical  method  is  to  be 
preferred  for  calculating  the  above.  In  the  past  the  preference  has 
been  for  the  trapezoidal  rule,  mainly  because  of  its  simplicity. 

Clements  (8)  not  only  used  the  trapezoidal  rule  but  also  Filon's  rule  in 
his  research  work  and  makes  an  interesting  comparison  of  the  two.  He 
favours  the  use  of  Filon’s  method  but  asserts  that  this  method  may 
require  inconveniently  long  computing  time. 

Recently,  a  highly  efficient  algorithm  for  calculating  the 
Fourier  transform  of  sampled  data  series  has  been  suggested  by  Cooley 
and  Tukey  (12).  Their  fast  Fourier  transform  (FFT)  technique  has  been 
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lauded  as  a  breakthrough  by  workers  in  the  field  of  digital  frequency 
analysis  (3) . 


A  discussion  of  each  of  the  three  methods,  experimental  work 
highlighting  the  differences  in  the  techniques,  and  a  concluding 
section  comparing  their  relative  worth  follows: 

Trapezoidal  Rule: 

This  is  the  simplest  method,  the  formula  employed  to  approxi¬ 
mate  the  integral  being 


fha  f (t)  dt  ~  At  [  £o  +  f-L  +  f2  + 


11-1 


and  is  seen  in  Figure  12  to  represent  the  summation  of  the  area  of  the 
trapezoids  inscribed  under  the  curve  y  =  f(t). 


y 


t 


Fig.  12  Trapezoidal  Approximation  of  a  Continuous  Function 
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By  simple  observation  it  is  apparent  that  the  less  "curved"  the 
function  y(t),  the  better  the  approximation.  Unfortunately,  the  quad¬ 
rature  products  to  be  evaluated  either  have  a  sin  mt  or  cos  mt  term,  and, 
as  a  result,  accuracy  decreases  as  frequency  increases.  For  a  further 
discussion  of  the  trapezoidal  rule  see  Appendix  II-A. 

Filon's  Method: 

Filon  (22) ,  interested  in  the  computation  of  trigonometric 
integrals,  developed  a  quadrature  method  which  circumvents  the  problem 
of  curve  fitting  highly  oscillating  functions  such  as  f(t)  sin  cot  and 
f(t)  cos  cot.  Basically,  the  time  function  is  approximated  by  a  quad¬ 
ratic  as  in  Simpson's  rule,  the  range  of  interest  being  divided  into 
an  even  number  of  increments;  then  the  integral  is  evaluated  analytically. 
The  method,  therefore,  does  not  break  down  due  to  high  frequencies. 
However,  the  method  still  entails  a  numerical  approximation  and  there¬ 
fore,  of  course,  is  still  subject  to  error. 


Consider  a  curve  divided  into  an  equal  number  of  increments,  At 


wide.  Denote  the  dividing  points  by  fQ,  f^, 
f(a)  sin  ma  and  f0  =  f(b)  sin  mb.  Then, 


f n  where  f „ 
2n  0 


/  b  f(t)  sin  cot  dt  ~  At  {  a[f(a)  cos  ma  -  f(b)  cos  mb] 
a 


11-2 
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where 

0  =  ui  At 

a  =  _1  +  sin  26  -  2  sin28 

8  20 2  03 

6  =  2  / cos26  +  1  -  sin  26  \ 

\  e2  e3  / 

y  =  4  /  sin  6  -  cos  6  \ 

\  63  62  / 


Also , 


f  f(t)  cos  oot  dt  "  At  (a[f(a)  cos  (wa  +  tt/2)  -  f(b)  cos  (cob  +  tt/2)  ] 
a 


+  gs2n  +  ts2  n  _  1  } 


...  11-3 


where  now 

fg  =  f(a)  cos  wa 

f„  =  f(b)  cos  u)b. 

2n 

It  should  be  mentioned  that,  when  0  tends  to  zero,  a,  6,  and  y  must  be 
calculated  from  series  expansion  formulae  to  prevent  loss  of  significant 
digits,  as  shown  in  Appendix  II-B.  The  obvious  disadvantage  of  this 
method  compared  with  the  trapezoidal  rule  is  the  extra  computation,  and 
possibly  extra  roundoff  error,  involved. 
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The  Fast  Fourier  Transform  Method: 

In  1965,  Cooley  and  Tukey  (12)  described  a  highly  efficient 
algorithm  for  computing  the  Fourier  transform  of  a  time  sampled  series. 
This  method,  now  known  as  the  fast  Fourier  transform  (FFT) ,  has  al¬ 
ready  proved  to  be  extremely  useful  in  a  number  of  applications  (2» 

11,  51) .  It  is  based  on  the  fact  that  the  calculation  of  the 
Fourier  coefficients  can  be  carried  out  iteratively,  resulting  in 
considerable  saving  of  computational  effort.  If,  for  example,  a  time 
series  consists  of  N  =  2n  samples,  then  the  number  of  arithmetic 
operations  required  to  evaluate  all  N  associated  coefficients  is 

approximately  2nN  =  2N  log2  N.  In  comparison,  the  number  of  opera- 

2 

tions  required  by  a  straightforward  procedure  is  N  .  Therefore,  the 
theoretical  speed  ratio  is 


S.R.  =  N_ 

2n 

Consider  a  series  for  which  N  =  1024  =  2  then  the  S.R.  ~  50;  this 
means  that  the  FFT  only  requires  about  2  per  cent  of  the  time  normally 
required  for  the  evaluation  of  the  Fourier  transform  of  a  1024-data 
point  series. 

Unlike  the  two  previous  quadrature  methods  already  discussed, 
where  the  frequency  values  may  be  arbitrarily  chosen  by  the  user,  the 
frequency  range  when  using  the  FFT  is  fixed  by  the  relation 

oj,  =  2tt i  ,  i  =  1,  2,  . ,N  ...11 

1  NAt 


■ 


'  ■  • 
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Thus  the  minimum  frequency  is 


co 

min 


2  TT 
NAt 


. ..  11-5 


and  the  maximum  frequency  is 


a)  =  2tt 

max  - — 

At 


. ..  11-6 


In  Chapter  10  it  was  pointed  out  that  aliasing  of  a  sampled  data 
series  occurs  at  the  Nyquist  frequency,  i.e.. 


TT 

At 


This  indicates  that  half  of  the  frequency  points  calculated  by 
the  FFT  suffer  from  the  effects  of  aliasing,  and  only  the  points  up  to 
N/2  should  be  regarded  as  valid.  The  limitation  of  the  frequency  range, 
and  the  fact  that  half  the  frequency  points  calculated  are  invalid,  are 
two  disadvantages  of  this  technique.  The  reader  is  referred  to 
Appendix  II- C  for  a  description  of  the  mathematical  basis  of  the  FFT 
algorithm. 


11.2  Comparison  of  Methods 

A  comparative  study  of  the  three  data  analysis  methods  was  under¬ 
taken.  Of  especial  interest  was  the  sensitivity  of  the  methods  to  loss 
of  data  points,  and  increasing  frequency. 
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FIG.  13  NORMALIZED  FREQUENCY  CONTENT  OF  RAMP  PULSE 
USED  IN  RUNS  13  AND  14 
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FIG.  14  COMPARISON  OF  AMPLITUDE  RATIOS  AS  GIVEN  BY  EACH 
OF  THE  DATA  ANALYSIS  METHODS  FOR  RUN  13 
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FIG.  15  COMPARISON  OF  PHASE  ANGLES  AS  GIVEN  BY  EACH 
OF  THE  DATA  ANALYSIS  METHODS  FOR  RUN  13 
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FIG.  16  COMPARISON  OF  AMPLITUDE  RATIOS  AS  GIVEN  BY  EACH 
OF  THE  DATA  ANALYSIS  METHODS  FOR  RUN  14 
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FIG.  17  COMPARISON  OF  PHASE  ANGLES  AS  GIVEN  BY  EACH 
.  OF  THE  DATA  ANALYSIS  METHODS  FOR  RUN  14 
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FIG.  18  NORMALIZED  FREQUENCY  CONTENT  OF  RAMP  PULSE 
USED  IN  RUNS  15  AND  16 
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FIG.  19  COMPARISON  OF  AMPLITUDE  RATIOS  AS  GIVEN  BY 
OF  THE  DATA  ANALYSIS  METHODS  FOR  RUN  15 
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FIG.  20  COMPARISON  OF  PHASE  ANGLES  AS  GIVEN  BY  EACH 
OF  THE  DATA  ANALYSIS  METHODS  FOR  RUN  15 
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FIG.  21  COMPARISON  OF  AMPLITUDE  RATIOS  AS  GIVEN  BY  EACH 
OF  THE  DATA  ANALYSIS  METHODS  FOR  RUN  16 
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FIG.  22  COMPARISON  OF  PHASE  ANGLES  AS  GIVEN  BY  EACH 
OF  THE  DATA  ANALYSIS  METHODS  FOR  RUN  16 
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frequency  points,  the  two  quadrature  methods  took  between 
3-1/2  and  4  minutes  to  execute.  In  comparison,  the  FFT 
took  only  1-1/2  minutes  to  calculate  all  of  its  512 
associated  frequency  points.  This  indicates  a  speed  ratio 
of  over  60. 

(ii)  Unlike  the  quadrature  method,  the  accuracy  of  the  FFT  is 
not  dependent  on  a  curve-fitting  technique,  and  as  such  is 
not  affected  by  the  higher  frequencies.  Due  to  aliasing, 
however,  at  frequencies  greater  than  w  =  n/At,  the  method 
is  inaccurate. 

(iii)  The  accuracy  of  the  method  suffers  no  more  than  the 
trapezoidal  rule  and  even  less  than  Filon’s  method,  as  the 
number  of  data  points  is  reduced. 

(iv)  For  the  system  conditions  chosen,  the  low  frequency 
amplitude  ratio  and  phase  angle  values  differ  noticeably 
from  their  respective  theoretical  curves.  For  example,  in 
Run  13,  the  amplitude  ratio  at  the  corner  frequency, 
instead  of  being  0.707,  is  approximately  0.9,  and  the 
phase  angle  result  is  similarly  in  error.  However,  the 
calculated  results  quickly  assume  the  theoretical  value 

as  the  Nyquist  frequency  is  approached.  Thereafter  the 
FFT  results  deviate  sharply  from  the  true  curve. 


•  V: 
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(b)  Filon’s  Method: 

(i)  Of  the  three  data  analysis  methods  this  is  the  most 
sensitive  to  loss  of  data  points,  as  shown  by  the  serious 
deviation  of  the  computed  frequency  results  from  the 
theoretical  curves  which  occur  in  Runs  14  and  16.  Run  16 
does  show,  however,  that  the  sensitivity  of  Filon’s  method 
to  loss  of  data  points  is  partly  offset  by  its  insensitivity 
to  high  frequencies. 

(ii)  Much  more  programming  effort  is  required  for  this  method 
as  compared  to  the  trapezoidal  rule  technique. 

(c)  Trapezoidal  Rule  Method: 

(i)  At  the  lower  frequency  ranges,  the  accuracy  of  this  method 
is  better  than  predicted  by  the  Lees  and  Dougherty  criteria. 
As  expected,  however,  the  results  deteriorate  as  frequency 
is  increased. 

(ii)  In  all  cases,  the  amplitude  ratio  values  calculated  at  the 
low  end  of  each  frequency  range  were  very  accurate,  allow¬ 
ing  an  excellent  estimation  to  be  made  of  the  steady  state 
gain  and  time  constant  of  the  system. 

A  limitation  of  the  above  study  is  that  the  effect  of  pulse  shape 
on  the  data  analysis  methods  is  not  known.  Even  so,  some  overriding 


conclusions  may  be  made,  as  follows: 


' 
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(a)  The  difficulty  of  writing  or  acquiring  an  FFT  program  may  well  be 
compensated  for  by  its  high  computational  efficiency  and  speed; 
care  must  be  taken,  however,  to  use  a  sufficiently  high  number  of 
time  data  points  and  small  value  of  sampling  interval,  to  ensure 
accurate  low  frequency  information. 

(b)  If  a  quadrature  method  is  used,  the  trapezoidal  rule  should  be 
satisfactory  for  most  engineering  work.  In  this  study  the 
trapezoidal  method  not  only  gives  equivalent,  but  in  most  cases 
better  results  than  Filon’s  method,  and  it  also  involves  far  less 
programming  effort.  At  frequency  ranges  higher  than  those 
studied  in  this  work  it  is  expected  that  the  trapezoidal  rule 
will  not  be  as  accurate  as  the  Filon  method. 

(c)  The  calculated  amplitude  ratios  for  either  quadrature  method 
were  accurate,  in  most  cases,  beyond  those  frequencies  predicted 
by  the  Lees  and  Dougherty  criteria.  It  is  felt,  however,  that 
this  criteria  should  still  be  used  to  specify  a  time  data 
sampling  rate,  especially  in  any  practical  application. 
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12  CONSIDERATIONS  OF  PULSE  SHAPE  AND  WIDTH 


Prior  to  an  actual  pulse  test  there  is  the  question  of  what  form 
of  forcing  function  should  be  used  in  order  to  maximize  valid  frequency 
response  information.  It  is  known  that,  for  a  given  pulse,  as  the 
magnitude  of  the  Fourier  transform  amplitudes  decreases,  the  more 
difficult  it  is  to  extract  reliable  frequency  information.  Therefore, 
one  main  consideration  in  selecting  a  forcing  function  is  the  behaviour 
of  its  frequency  content  characteristic.  Previous  workers,  notably 
Lees  (36),  Dreifke  (17)  and  Clements  (8),  have  made  a  comprehensive 
study  of  what  constitutes  a  favourable  forcing  function  as  regards  to 
shape,  width  and  height,  and  the  following  conclusions  and  comments  are 
a  result  of  their  findings.  Examples  of  regular  pulse  shapes  are 
shown  in  Figure  23. 


Rectangular 


Triangular 


Ramp 


Rectangular  Doublet 


Fig.  23  Examples  of  Regular  Pulse  Shapes 
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To  facilitate  the  discussion  it  is  convenient  to  classify  the 
pulses  as  either  one-sided,  like  the  triangular,  or  double-sided,  like 
the  rectangular  doublet.  The  following  considerations  pertain  to  one¬ 
sided  type  pulses. 

A  wide  pulse,  with  a  duration  10  times  the  dominant  time  constant 
of  the  process,  although  giving  good  steady  state  gain  information, 
gives  relatively  poor  high  frequency  information.  As  the  pulse  width 
is  decreased,  the  frequency  content  of  the  forcing  function  increases, 
with  a  resulting  improvement  in  the  accuracy  of  the  computed  results  at 
the  higher  frequencies.  If  the  pulse  width  is  decreased  to  the  limit, 
a  pure  Dirac  impulse  will  be  realized.  This  function  exhibits  a  con¬ 
stant  frequency  spectrum  and  is  therefore  the  ideal  pulse  shape. 
Unfortunately,  a  forcing  function  of  infinitesimal  width  is  not  attain¬ 
able  experimentally.  Another  practical  aspect  is  that  too  narrow  a 
pulse  produces  such  an  insignificant  change  in  the  system  output  that 

the  latter  is  indiscernible  from  process  noise.  The  smallest  pulse 
width  generally  found  compatible  with  experimental  limitations  is 
approximately  one-tenth  of  the  major  time  constant  of  the  system.  On 
the  basis  of  the  above  considerations,  a  safe  guideline  when  initiat¬ 
ing  plant  tests  is  to  make  the  pulse  width  of  the  order  of  the  estimated 
dominant  time  constant  of  the  process.  The  greater  the  magnitude  of 
the  pulse  for  a  given  width,  the  higher  the  energy  transfer  into  or 
out  of  a  process.  Care  must  be  taken  to  ensure  that  the  pulse  is  not 
so  large  as  to  saturate  the  system,  but  not  too  small  as  to  produce  no 
significant  change  in  the  process  response. 


' 
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As  already  indicated,  the  frequency  content  characteristics  of  a 
pulse  are  profoundly  affected  by  its  shape.  In  general,  one-sided 
pulses  exhibit  their  maximum  frequency  content  at  zero  frequency,  and 
those  which  are  symmetrical  exhibit  one  or  more  zeros.  The  low  extreme 
of  frequency  content  is  exhibited  by  the  rectangular  pulse,  the  high 
extreme  by  the  impulse  function.  These  factors  are  reflected  in  the 
increasing  frequency  content  of  pulses  which  move  away  from  the  rect¬ 
angular,  f lat-in-the-middle  shape,  to  forms  exhibiting  a  sharp  central 
spike.  Two-sided  equal-weighted  pulses  such  as  the  rectangular  doublet 
have  a  frequency  content  of  zero  at  zero  frequency  and  exhibit  a  maxi¬ 
mum  at  some  higher  frequency.  This  implies  that  the  double-sided  pulse 
may  be  very  useful  in  maximizing  system  excitation  in  selected 
frequency  regions. 

In  an  actual  plant  test,  the  implementation  of  a  regular  shape 
forcing  function  will  usually  not  be  practical,  the  pulse  probably 
being  introduced  by  means  of  a  hand  valve.  However,  when  analyzing 
the  data,  if  the  frequency  content  of  the  input  function  is  calculated, 
regions  of  uncertainty  in  the  frequency  diagram,  due  to  the  Fourier 
amplitudes  approaching  zero,  may  be  assessed.  The  frequency  results 
are  usually  assumed  to  be  invalid  when  the  normalized  frequency  drops 
below  0.25  to  0.3  (8).  In  summary,  the  pulse  chosen  for  a  particular 
test  should  be  a  shape  and  width  as  to  ensure  that  it  contains  a  wide 
enough  band  of  component  frequencies  so  that  the  system  dynamics  are 
excited  over  the  frequency  range  of  interest. 
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13  GUIDELINES  TO  DESIGNING  A  PULSE  TESTING  EXPERIMENT 


On  the  basis  of  previous  theoretical  studies,  it  is  possible  to 
outline  a  procedure  for  designing  a  pulse  test.  It  is  hoped  that  the 
following  guidelines  will  be  useful  to  anyone  wishing  to  make  a  dynamic 
analysis  of  a  process  or  piece  of  equipment.  Prior  to  any  design, 
however,  there  are  questions  concerning  the  actual  experimental  equip¬ 
ment  which  must  be  satisfactorily  answered  before  a  pulse  test  can  be 
expected  to  reveal  accurate  frequency  information: 

(a)  Are  the  sensing  devices  accurate,  sufficiently  sensitive,  and 
actually  measuring  the  desired  signals? 

(b)  Can  a  pulse  be  implemented  safely  and  without  deleterious  effects 
on  the  end  products  of  the  process? 

(c)  A  closed  pulse  is  assumed  in  the  theoretical  development  of  the 
data  analysis  methods.  Can  the  forcing  function  be  brought  back 
within  a  reasonable  proximity  of  the  steady  state  condition 
existing  prior  to  the  pulse? 

(d)  Are  the  recording  devices  satisfactory?  For  example,  if  ordinary 
pen  recorders  are  to  be  used,  is  the  chart  speed  high  enough,  is 
best  use  being  made  of  the  span  of  the  instrument,  and  is  it 
zeroed  out  and  calibrated  properly?  If  an  on-line  system  is  used, 
are  the  analog  inputs  successfully  shielded  from  extraneous 
noise,  are  the  conversion  programs  correct,  and  is  the  data  log¬ 
ging  being  done  at  the  proper  sampling  intervals? 


C.  r.  ,  ,  L1  J  ..Xiir  ..  »  ] 
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If  these  questions  are  answered  satisfactorily,  consideration  may 
be  given  to  the  design  of  the  pulse  test.  It  is  suggested  that  the  steps 
in  the  design  procedure  be  followed  in  the  order  given  below.  It  is 
assumed  that  the  experimenter  has  no  prior  knowledge  at  all  of  the 
dynamics  of  the  system  to  be  studied. 

(i)  For  convenience,  assume  that  the  process  under  consideration 
exhibits  the  behaviour  of  a  first-order  system.  For  most  chemical 
processes  this  is  not  too  unreal  an  assumption. 

(ii)  Estimate  the  time  constant  for  this  supposedly  first-order  process. 
This  may  be  done  in  a  number  of  ways .  Previous  operating  records 
will  probably  furnish  information  on  step  changes  in  the  process. 
Dividing  the  length  of  time  the  response  variable  took  to  reach  a 
new  steady  state  value  by  5  will  give  an  approximate  value  for  t. 

If  a  continuous  flow  system,  another  simple  method  is  to  divide 
the  hold-up  of  the  equipment  between  the  point  where  the  pulse  is 
applied  to  where  the  response  is  measured,  by  the  volumetric  flow 
rate . 

(iii)  The  intended  pulse  should  have  a  width  equal  to  x.  Also,  using 
the  knowledge  of  previous  operating  conditions,  it  should  be  of 
sufficient  magnitude  to  cause  a  significant  change  in  the  response 
without  saturating  the  system.  Finally,  a  sharp  pointed  pulse, 
such  as  the  triangular  or  ramp  function,  should  be  implemented, 

if  possible,  due  to  their  favourable  frequency  content  charac¬ 
teristic.  This,  however,  is  not  essential  for  engineering  design 
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purposes,  as  even  the  rectangular  pulse  will  usually  provide 
sufficiently  accurate  frequency  response  information. 

(iv)  Prior  to  testing,  the  experimenter  will  be  unsure  of  the 

frequency  range  over  which  he  should  calculate  the  amplitude 

ratio  and  phase  angle.  In  the  absence  of  known  upper  and  lower 

frequency  limits  it  is  suggested  that  the  initial  range  of  interest 

be  taken  as  one  decade  above  and  below  the  corner  frequency  (wc) » 

a)  =  1/t. 

c 

(v)  To  avoid  truncation  errors,  the  time  histories  should  be  recorded 
until  the  system  returns  to  the  steady  state  condition.  For  con¬ 
venience,  it  may  be  assumed  that  this  is  reached  after  a  period  of 
five  time  constants  after  the  completion  of  the  pulse.  Allowing 

a  period  of  x  prior  to  the  pulse,  to  ensure  the  accurate  recording 
of  the  steady  state  values,  the  total  time  record  for  the  whole 
pulse  test  should  therefore  be  equal  to  approximately  seven  time 
constants . 

(vi)  Using  the  maximum  value  in  the  frequency  range  of  interest,  and 
the  Lees  and  Dougherty  criteria  for  avoiding  numerical  approxima¬ 
tion  errors,  the  sampling  interval  may  be  computed  from 

At  =  0. In 

CO 

max 

Dividing  this  value  into  the  required  length  of  the  time  records 
yields  the  number  of  data  points . 
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With  the  above  information  a  pulse  test  may  be  designed.  Upon 
analyzing  the  resulting  data,  the  frequency  content  of  the  input  pulse 
should  be  calculated.  This  aids  in  determining  the  validity  of  the 
frequency  response  diagram.  Knowledge  gained  from  the  first  test  should 
give  an  estimate  of  the  order  of  the  system  studied  and  the  value  of  the 
dominant  time  constant.  If  further  tests  are  deemed  necessary  they  can 
be  readily  designed  on  the  basis  of  the  information  obtained  from  the 
initial  test. 

It  should  be  mentioned  that,  on  analyzing  the  time  data,  it  is 
possible  that  "scatter"  and  therefore  doubtfully  valid  frequency  results 
may  occur  at  a  frequency  content  greater  than  the  usually  cited  value 
of  0.3.  In  this  case  the  experimenter  should  not  dismiss  the  possibility 
that  the  erroneous  results  may  not  only  be  due  to  one  or  more  of  the 
causes  of  inaccuracies  presented  earlier,  but  may  also  be  due  to  the 
presence  of  extraneous  input  disturbances  to  the  process  during  the 
pulse  test.  If  this  is  so,  measures  may  be  warranted  to  eliminate,  or 
if  this  is  not  possible,  to  reduce  the  magnitude  of  these  input 


disturbances . 
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PART  III 


ON-LINE  PULSE  TESTING 


i 

. 


14  DETAILS  OF  EXPERIMENTAL  INVESTIGATION 


14.1  Introduction 

The  experimental  work  was  carried  out  on  a  double  pipe,  coun¬ 
tercurrent,  liquid-liquid  heat  exchanger,  operated  in  a  DDC  mode 
using  an  IBM  1800  computer.  This  investigation  continues  the  pulse 
testing  work  of  R.  S.  Lees  (35),  who  designed  the  experimental 
apparatus,  did  the  necessary  computer  interfacing,  wrote  the  loop 
records  required  for  the  DDC  work  and  the  programs  necessary  for  the 
shell  thermopile  temperature  and  tube  fluid  flow  conversion. 

The  investigation  may  be  divided  into  three  main  areas:  the 
equipment  and  its  operation,  the  process  programs  to  initiate  the 
desired  pulse  and  to  retrieve  and  store  the  time  data,  and  the  non¬ 
process  programs  to  perform  the  data  analysis  and  plot  the  results. 
In  this  chapter  the  heat  exchanger  and  ancilliary  equipment  required 
for  the  on-line  pulse  testing  are  described;  details  of  the  process 
and  nonprocess  programs  are  given  in  Chapter  15. 


14.2  Description  of  Experimental  Equipment 

Figure  24  shows  the  location  of  the  sensing  and  control  devices, 
the  arrangement  of  the  analog  input-output  points,  and  the  general 
layout  of  the  test  heat  exchanger.  Figure  25  shows  the  information 
flow  scheme  between  the  apparatus  and  the  computer  necessary  for 
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FIG.  24  LOCATION  OF  ANALOG  INPUT-OUTPUT  POINTS  AND  GENERAL  LAYOUT 
OF  EXPERIMENTAL  HEAT  EXCHANGER  AND  AUXILIARY  EQUIPMENT 
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direct  digital  control  and  data  acquisition. 


The  test  heat  exchanger  was  constructed  of  two  concentric 
copper  tubes  10  feet  in  length.  The  inside  tube  had  an  inside 
diameter  of  0.815  inches  and  a  wall  thickness  of  0.030  inches*  The 
inside  diameter  of  the  outside  tube  (shell)  was  1.541  inches  and  it 
had  a  wall  thickness  of  0.042  inches.  Spacers  placed  along  the 
heat  exchanger  ensured  the  concentricity  of  the  tube  and  shell  as 
well  as  preventing  sagging  in  the  unsupported  sections. 

The  IBM  1800  computer  is  a  32-K  word  machine  with  a  fixed 
length  word  of  16  bits.  Three  disks,  each  of  512-K  word  capacity, 
provide  the  user  with  random  access  storage.  Peripheral  devices  for 
input-output  operations  include  an  IBM  1442  card-read-punch  machine, 
two  IBM  1816  keyboard-printer  typewriters,  an  IBM  1132  printer,  an 
IBM  1627  plotter,  and  a  Tektronix  Incorporated  Type  611  oscilloscope. 

Although  not  indicated  in  Figure  24,  the  temperature  of  the  hot 
water  fed  to  the  shell  side  of  the  test  heat  exchanger  was  also 
a  controlled  variable.  The  hot  process  water  was  maintained  at  the 
desired  temperature  by  means  of  a  heat  exchanger  (two-tube  pass,  one- 
shell  pass).  The  flow  of  the  shell  side  fluid  was  controlled  using 
a  3/4-inch  Fisher  Governor  control  valve  with  linear  trim.  A  similar 
valve  was  used  to  control  the  cold  process  water  in  the  tube  of  the 
heat  exchanger.  Upstream  liquid  pressures  were  maintained  by  use  of 
Weinman  centrifugal  pumps,  and  liquid  flows  were  measured  by  turbine 
flowmeters  manufactured  by  Flow  Technology  Incorporated.  The 
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measurable  range  of  flow  on  the  tube  side  was  1.5  to  22.0  U.S.  gallons 
per  minute,  and  on  the  shell  side  1.5  to  40.0  U.S.  gallons  per  minute. 

For  the  purpose  of  heat  balance  calculations,  four  miniature 
Thermoelectric  iron/constantan  thermocouples  were  located  at  the  two 
stream  inlets  and  outlets  of  the  heat  exchanger.  Another  thermo¬ 
couple,  chromel/constantan,  was  situated  in  an  oil  bath  which  acted 
as  a  cold  junction  reference  for  the  thermopile.  The  latter  consisted 
of  five  copper/constantan  thermocouples  in  series  and  was  used  to 
measure  the  response  variable,  namely  the  shell  fluid  outlet  tempera¬ 
ture.  The  forcing  variable,  tube  liquid  flow,  was  measured  by  the 
pertinent  turbine  flowmeter. 

To  achieve  data  acquisition  and  direct  digital  control  of  a 
process,  information  must  flow  between  the  computer  and  the  sensing 

and  control  elements.  Figure  25  shows  the  information  flow  scheme 

/ 

pertinent  to  this  work.  The  analog  signals  from  the  sensing  devices, 
i.e.,  thermopile,  thermocouples,  and  turbine  flowmeters,  are  fed  to 
a  bank  of  multiplexer  points.  Once  every  second  the  computer  scans 
a  point  by  closing  a  contact,  thereby  allowing  the  signal  to  be  read 
by  the  analog-to-digital  converter  (ADC).  Here  it  is  converted  to 
an  integer  value  acceptable  to  the  computer,  ranging  from  0  to  32767. 
Signals  going  from  the  computer  to  the  process  have  to  be  converted 

to  an  analog  form;  this  is  done  by  the  digital-to-analog  converter 

(DAC).  A  signal  emanating  from  the  DAC  is  directed  to  its  proper 

current  output  station  by  means  of  closing  a  switch  (a  pulse  output 

point).  The  current  output  station  is  a  zero  order  hold  amplifier 
with  direct  access  to  a  final  control  element.  The  output  signal  to 
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the  control  device  is  held  constant  by  the  current  output  station 
until  a  new  signal  is  received  from  the  DAC.  The  scanning  and 
switching  operation  for  both  analog  inputs  and  outputs  is  controlled 
automatically  by  supervising  software. 

A  more  detailed  description  of  the  equipment  and  its  operation 
is  given  by  Lees  (35). 


' 


100 


15  COMPUTER  PROGRAMS 


15.1  Introduction 

For  the  experimental  and  theoretical  analysis  work  advantage 
was  taken  of  the  real  time  capabilities  of  the  IBM  1800  computer. 
Programs  were  written  to  attain  the  following  objectives: 

(a)  Implementation  of  on-line  pulse  testing  of  the  heat  exchanger 

(b)  Acquisition,  conversion,  storage  and  display  of  the  resulting 
time  data 

(c)  Processing  of  the  data  with  the  options  of  data  point  reduc¬ 
tion  and  method  of  analysis 

(d)  Output  of  the  frequency  results  in  a  convenient  format. 

Two  sets  of  linked  programs  were  written,  one  for  those  con¬ 
cerned  with  the  process,  (a)  and  (b)  above,  and  one  for  those  con¬ 
cerned  with  the  time  data  analysis,  (c)  and  (d)  above.  Each  program, 
either  process  or  nonprocess,  had  to  be  small  enough  to  fit  into  the 
variable  core  area  of  the  computer;  in  this  case  9.7-K  words.  This 
placed  a  heavy  emphasis  on  efficient  disk  usage  and  program  linking, 
and  accounts  for  the  style  of  program  writing.  A  detailed  descrip¬ 
tion  of  the  two  sets  of  programs  follows. 


. 
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15.2  Process  Programs 

These  programs  allow  the  user  to  make  a  pulse  test  on  the  heat 
exchanger,  concurrent  with  high  speed  data  acquisition  of  the  tube 
flow  (forcing  variable)  and  shell  fluid  outlet  temperature  (response 
variable).  Upon  completion  of  the  pulse  test,  the  time  data,  being 
in  integer  format,  are  converted  to  engineering  values,  which  are 
then  placed  in  user-specified  disk  files.  Output  of  the  data  is 
possible  on  the  oscilloscope,  the  X-Y  plotter,  or  on  punched  cards. 

A  flowchart  of  the  programs  is  given  in  Figure  26  ,  and  a 
general  description  of  the  function  of  each  of  the  process  programs 
follows.  A  listing  of  the  programs  appears  in  Appendix  III-A. 

HTBAL : 

This  program  does  an  on-line  heat  balance  of  the  heat 
exchanger  over  a  given  length  of  time.  The  data  acquisition  period 
is  entered  by  the  user  via  a  keyboard  entry  request  feature.  Every 
scan  time  (1  second)  the  following  variables  are  recorded:  tube  flow, 
shell  flow, and  the  temperatures  at  the  tube  inlet,  tube  outlet,  shell 
inlet  and  shell  outlet.  The  values  are  averaged  and  a  heat  balance 
calculated,  the  results  being  printed  out  on  the  typewriter.  On  the 
basis  of  this  information  the  user  may  abort  the  run  (data  switch  5) 
or  call  in  the  next  program  (data  switch  6). 


GENER: 

GENER  is  a  pulse  generation  program;  it  consists  of  a  mainline 
and  several  servicing  subroutines.  Each  subroutine  generates  a 
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different  pulse  shape,  i.e.. 

Subroutine  1  -  Rectangular  pulse 
Subroutine  2  -  Ramp  pulse 
Subroutine  3  -  Triangular  pulse 
Subroutine  4  -  Displaced  cosine  pulse 
Subroutine  5  -  Rectangular  doublet  pulse 

The  mainline  requests  the  user,  via  the  keyboard,  to  specify  the  fol¬ 
lowing  parameters: 

ISUB  -  The  number  of  the  subroutine  which  generates  the 
desired  pulse  shape 
IPULW  -  The  pulse  width  in  seconds 

IRATE  -  The  number  of  times  per  second  the  control  valve  will 
be  pulsed. 

* 

The  pulse  values,  Y^,  are  then  generated  between  zero  and  plus  or 
minus  one.  In  the  actual  pulsing  program,  these  values  are  converted 
to  controller  output  values  (%)  acceptable  to  the  control  valve.  It 
should  be  noted  that,  in  all  the  subroutines  for  pulse  generation, 
five  seconds  of  steady  state  values  precede  the  pulse. 


TPLOT : 

To  ensure  that  the  subroutine  called  is  the  desired  one  and 
that  the  subroutine  does  actually  generate  the  correct-shaped  forcing 
function,  this  program  automatically  displays  the  generated  pulse  on 
the  oscilloscope.  At  this  stage  the  user  may  request  the  pulse  values 
to  be  plotted  on  the  X-Y  plotter  so  that  a  permanent  record  may  be 


obtained . 


’ 


. 


' 
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PULSE: 


This  program  has  two  main  functions:  to  pulse  the  control 
valve  using  the  generated  values  from  GENER,  while  concurrently 
performing  high  speed  data  acquisition  of  the  forcing  and  response 
variables.  The  program  first  requests  the  span  of  the  pulse  in  per 
cent,  and  this  value  is  entered  via  the  keyboard.  The  pulse  values 
generated  in  GENER  are  then  modified  using  the  relationship: 


Y  x  SPAN  (%) 


present  output  to 
control  valve  (%) 


The  flow  control  loop  is  placed  on  "manual",  then  the 
implementation  of  the  pulse  is  begun,  the  signals  being  sent  to  the 
valve  at  a  rate  of  5  times  per  second;  at  the  same  time  the  high 
speed  data  acquisition  subroutines  are  activated.  These  routines 
collect  data  at  40  points  per  second  for  40  seconds.  Each  set  of 
1600  values  is  stored  on  disk.  On  completion  of  the  40-second  run, 
of  which  the  first  5  seconds  are  steady  state  data,  the  flow  con¬ 
trol  loop  is  put  back  on  "automatic". 


ASORT: 


Employing  the  system  subroutine  FLSRT,  this  program  takes 
the  data  collected  during  the  high  speed  data  acquisition  period  and 
sorts  it  into  user-specified  files. 


VTFL0: 


The  flow  data  is  in  integer  form,  each  value  being  between  0 


and  32767.  To  be  meaningful  this  raw  data  has  to  be  converted  to 
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engineering  units.  The  VTFLO  program  converts  the  raw  integer  flow 
data  into  values  of  U.S.  gallons  per  minute.  Each  raw  data  point  is 
converted  to  an  equivalent  voltage  value  (5-volt  input  point  at  the 
ADC).  Then,  using  a  conversion  equation  for  the  flow  turbine,  the 
voltage  value  is  converted  to  an  equivalent  flow  in  U.S.  gallons  per 
minute.  After  conversion,  the  engineering  data  is  stored  on  disk. 

MVTMP: 

This  is  the  conversion  program  for  the  temperature  data.  Be¬ 
fore  the  actual  conversion  can  take  place,  however,  cold  junction 
compensation  must  be  applied.  It  should  be  recalled  here  that  the 
thermopile  uses  an  oil  bath  for  its  cold  junction.  The  temperature 
of  the  oil  bath  is  given  by  a  chromel/constantan  thermocouple  which, 
in  turn,  has  for  its  cold  junction  the  reference  bulb  temperature 
of  the  computer.  The  reference  bulb  temperature  is  calculated  first. 
This  is  done  by  supplying  the  computer  reference  and  resistance 
bulb  voltages  to  a  given  system  equation.  Cold  junction  compensa¬ 
tion  is  now  made  to  the  chromel/constantan  thermocouple  situated  in 
the  oil  bath.  Knowing  the  correct  temperature  of  the  oil  bath,  cold 
junction  compensation  is  now  possible  for  the  outlet  shell  tempera¬ 
ture  thermopile.  The  raw  integer  data  from  the  thermopile  is  con¬ 
verted  to  equivalent  millivoltage  values  (20  mv  input  at  ADC)  using 
a  copper /constantan  conversion  equation.  These  values  are  converted 
to  temperature  data  (degrees  Fahrenheit)  and  stored  in  disk  files. 


BPLOT : 

This  program  allows  the  user  to  plot  out  the  two  records  of 
time  data  on  the  oscilloscope,  thus  providing  immediate  visual 
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results  of  the  experimental  run.  If  a  permanent  record  is  desired, 
the  data  may  be  plotted  on  the  X-Y  plotter.  The  user  may  also 
specify  the  data  to  be  punched  out  on  cards  in  a  format  suitable 
for  the  nonprocess  analysis  programs. 


15.3  Nonprocess  Programs 

For  the  theoretical  studies  as  well  as  the  analysis  of  the 
experimental  time  data,  it  was  necessary  to  design  the  nonprocess 
programs  to  allow  certain  user-specified  options.  These  were: 

(a)  Data  point  reduction 

(b)  Data  analysis  method 

(c)  Acceptance  of  user-generated  frequency  results.  In  the 
theoretical  studies  this  allowed  the  true  frequency  response 
curve  of  the  process  to  be  plotted  on  a  Bode  diagram,  thereby 
facilitating  the  comparison  of  the  data  analysis  methods 

(d)  Printing  or  plotting  of  the  frequency  results 

(e)  Allowing  the  user  to  return  to  the  initial  program  for  a 

re-execution  of  the  series. 

To  do  this,  nine  linked  mainline  programs  were  written.  A 
flowchart  indicating  the  structure  of  the  programs  appears  in 
Figure  27.  A  detailed  description  of  each  program  follows,  and  a 
listing  appears  in  Appendix  III-B. 
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PTAP: 


This  is  the  executive  program.  It  permits  the  user,  either 
via  the  card  reader  or  the  typewriter,  to  enter  the  following 
parameters  and  flags: 


(a)  Data  type  and  source: 


(i) 

Time  data  -  on  cards 

) 

) 

) 

1 

(ii) 

Time  data  -  in  disk  files 

2 

) 

ITYPE 

= 

(iii) 

Frequency  results  -  on  cards 

) 

) 

) 

3 

(iv) 

Frequency  results  -  in  disk  files 

4 

(b) 

Desired 

output  for  the  frequency  results: 

(i) 

Write  out  only 

) 

) 

) 

) 

) 

1 

(ii) 

Plot  out  only 

I  OUT 

=  2 

(iii) 

Write  out  and  plot  out 

3 

(c) 

Device  for  data  parameter  entry: 

(i) 

Via  the  keyboard 

) 

1 

) 

IDEV 

25 

(ii) 

Via  the  card  reader 

) 

2 

(d) 

Is  another  program  execution  required? 

(i) 

Yes,  return  to  PTAP 

) 

1 

) 

IEXIT 

= 

(ii) 

No,  exit 

) 

2 

(e) 

Frequency  parameters: 

(i) 

WI  -  initial  frequency 

(ii) 

WF  -  final  frequency 

(iii) 

NNW  -  desired  number  of  frequency 

points 

to  be 

calcul 

ated  between  WI  and  WF 


- 


■ 
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(f)  Time  data  parameters: 


(i)  NN  -  total  number  of  time  data  points 

(ii)  ISS  -  number  of  steady  state  data  points 

(iii)  IPUL  -  number  of  time  increments  constituting  the  width 
of  the  pulse 

(iv)  DT  -  time  increment 


(g)  Quadrature  method  to  be  used: 

(i)  Trapezoidal  rule 

(ii)  Filon’s  method 

(iii)  Fast  Fourier  transform 


)  1 

) 

)  IQUAD  =  2 

> 

)  3 


(h)  Is  time  data  reduction  required? 

(i)  All  time  data  points  are  to  be  used  ) 

) 

(ii)  Every  2nd  data  point  is  to  be  used  ) 

.  I  RED 


1 

2 


(n)  Every  nth  data  point  is  to  be  used  ) 


n 


When  entering  the  above  parameters  and  flags,  the  following  should 
be  noted: 

(a)  In  specifying  the  frequency  range  of  interest,  WI  and  WF  should 
be  whole  decades,  i.e.,  0.001,  1.0,  100,  etc. 


(b)  The  maximum  number  of  frequency  points  that  may  be  specified  is 


A0,  i.e., 


(NNW) 


max 


40 
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(c)  If  the  parameters  and  flags  are  to  be  read  via  the  card  reader, 
the  following  format  applies: 

WI,  WF,  NNW  -  (2F10.4,  13) 

NN,  ISS,  IPUL,  DT,  IQUAD,  IRED  -  (314,  F7.3,  212) 


WMAIN: 

If  the  frequency  results  are  already  available,  either  in 

files  or  on  cards,  PTAP  links  to  WMAIN.  The  latter  has  two  functions: 

(a)  If  the  frequency  results  are  already  in  the  files,  the  program 
links  to  the  mainline  APEN. 

(b)  If  the  frequency  results  are  on  cards,  the  program  reads  the  data 
into  core,  then  writes  it  out  in  the  appropriate  files.  The 
mainline  APEN  is  then  called. 

The  file  specifications  are  as  follows: 


51 

-  W 

52 

-  S(W)n 

53 

-  A.R. 

(W) 

54 

-  PHI. 

(W) 

55 

-  x(t) 

56 

-  Y(t) 

PMAIN : 

PMAIN  is  an  executive-type  program  which  reads  the  time  data 
either  off  cards  or  the  disk  and  performs  certain  necessary  functions 
prior  to  the  Fourier  transform  calculations.  The  functions  performed 


are : 


.  '  '  < 

■ 


■ 
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(a)  Writing  the  time  data  into  the  appropriate  records  on  disk  (PREAD). 

(b)  Writing  out  the  specified  parameters  (PWRIT) . 

(c)  Time  data  reduction  (REDN) . 

(d)  Calculation  of  a  vector  of  frequency  points  (WVEC) . 

(e)  Calculation  of  the  steady  state  values  of  the  raw  input  data 
(STEDY) . 

(f)  Calculation  of  the  area  under  the  input  pulse  (STEDY). 

The  subroutines  which  service  the  mainline  are  shown  in  brackets. 
Depending  on  the  flag  IQUAD,  a  link  is  then  made  to  either  of  the 
three  quadrature  programs. 

TRAP: 

The  real  and  imaginary  parts  of  the  Fourier  transform  of  both 
the  input  and  the  output  time  data  are  calculated  using  the 
trapezoidal  rule.  The  resulting  vectors  are  then  written  on  the  disk. 


FILON : 

As  above,  except  that  Filon’s  method  is  used  to  calculate  the 
quadrature  products. 

FFT: 

As  above,  except  that  the  fast  Fourier  transform  method  is 
used.  The  program  used  here  was  supplied  by  Trinity  University, 
Texas . 


. 
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RSLTS : 

This  program  reads  the  quadrature  products  off  the  disk  and 
calculates  the  amplitude  ratio,  the  phase  lag  and  the  frequency 
content  at  each  desired  frequency  point. 

APEN: 

If  desired,  the  results  can  be  written  out  on  the  1443  printer, 
using  this  program.  Depending  on  the  flags  specified  in  PTAP,  one  of 
the  following  actions  is  taken  once  the  printing  has  been  executed. 

(a)  Link  to  the  plotting  program,  APLOT 

(b)  No  plot,  but  link  back  to  PTAP  for  another  execution 

(c)  Exit. 

APLOT : 

Plotting  of  the  results  is  achieved  by  APLOT.  Use  is  made  of 


the  plotting  subroutine  NWPLT.  On  completion  of  the  plot,  the  program 
either  links  back  to  PTAP,  or  exits. 


' 
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Fig.  26  FLOWCHART  SHOWING  STRUCTURE  OF  PROCESS  PROGRAMS 
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Fig.  26  Continued.... 
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VTFLO 

Conversion 
program  for 
flow  data 

1 

' 

MVTMP 


Conversion 
program  for 
temperature 
data 
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Output  of  time 
data  -  scope, 
plotter,  and 
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FIG  27 


FLOWCHART  OF  NONPROCESS  PROGRAMS 
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Fig.  27  Continued 


16  CHARACTERIZATION  OF  THE  TEST  HEAT  EXCHANGER 


16.1  Introduction 

The  double  pipe,  countercurrent,  liquid-liquid  heat  exchanger 
was  pulse  tested  on  line,  using  five  different  forcing  functions. 

The  pulses  employed  were: 

(a)  Rectangular 

(b)  Ramp 

(c)  Triangular 

(d)  Displaced  cosine 

(e)  Rectangular  doublet 

Symmetrical,  nonsymmetrical ,  one-sided  and  double-sided  pulses 
are  represented  in  the  above  group. 

The  time  data  obtained  during  the  dynamic  simulation  studies 
indicated  that  the  transient  response  of  the  shell  fluid  outlet 
temperature  to  the  tube  flow  forcing  could  be  represented  by  a  trans¬ 
portation  lag  and  first-order  system  in  series,  with  a  lag  of  1.0  and 
a  time  constant  of  2.5  seconds  respectively.  A  pulse  width  of  4 
seconds  duration  was  chosen  for  the  experimental  investigation.  This 
is  in  accordance  with  the  previously  made  recommendation  that  the 
width  of  a  forcing  function  should  be  of  the  same  order  of  magnitude 
as  the  dominant  time  constant  of  the  process. 

At  the  commencement  of  each  run,  the  heat  exchanger,  under 
direct  digital  control,  was  allowed  to  reach  a  steady  state  condition. 


■ 


■ 


■ 


Then,  by  means  of  a  keyboard  request  feature,  a  pulse  of  a  given 
shape,  width  and  magnitude,  was  specified  by  the  user.  It  was  pos¬ 
sible  to  change  the  output  signal  to  the  tube  flow  control  valve  at 
a  rate  of  five  times  per  second.  This  meant  that  a  flow  pulse  of 
four  seconds  duration  could  be  described  by  twenty  output  signal 
values.  This  was  found  to  be  a  sufficient  number  to  satisfactorily 
represent  the  pulse  shapes  employed  in  this  study. 

On  implementation  of  a  pulse,  the  two  time  histories  of 
interest,  the  tube  fluid  flow,  and  the  shell  fluid  outlet  tempera¬ 
ture,  were  recorded  as  two  sets  of  digitized  values.  Sixteen  hundred 
values  representing  a  forty-second  run  at  forty  points  per  second 
were  placed  in  disk  files  for  each  of  the  two  time  records.  Process 
programs  were  automatically  linked  to  convert  the  data  from  integer 
to  engineering  values.  For  presentation  purposes,  the  user,  by  means 
of  the  final  process  program,  BPLOT,  was  able  to  plot  out  the  time 
data  results  on  the  DACS  Centre’s  oscilloscope.  If  the  data  were 
desired  for  later  analysis,  the  results  were  plotted  on  the  1627 
plotter  and  also  punched  out  on  cards  in  a  format  acceptable  to  the 
nonprocess  data  analysis  program. 

Using  the  nonprocess  program,  the  normalized  frequency  con¬ 
tent,  amplitude  ratio  and  phase  angle  for  each  pulse  test  were 
calculated  over  the  frequency  range  0.01  to  10.0  radians  per  second. 
From  knowledge  of  the  frequency  content  characteristics  of  the 
pulses  used  in  the  study,  it  was  expected  that  accurate  frequency 
data  would  only  be  realized  up  to  approximately  3  radians  per  second. 
Using  this  figure,  the  Lees  and  Dougherty  (37)  criteria  indicated 


1 


r  • 
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the  necessity  to  use  only  20  data  values  per  second.  This  halved 
the  data  analysis  program  execution  time.  Computation  time  was 
also  significantly  reduced  by  using  only  the  initial  20  seconds 
of  each  time  record,  since  the  transients  had  usually  decayed  by 
the  end  of  this  period. 

The  trapezoidal  quadrature  method,  instead  of  the  FFT,  was 
used  in  the  data  analysis  because  of  steady  state  gain  considera¬ 
tions.  It  will  be  recalled  that  the  FFT,  although  a  more  efficient 
method,  does  not  allow  the  user  to  define  an  arbitrary  frequency 
range.  To  obtain  good  steady  state  gain  information,  the  low 
frequency  limit  of  0.01  radians  per  second  was  desired;  using  the 
FFT  the  lowest  frequency  possible  was  only  0.157  radians  per  second. 

The  experimental  conditions,  time  data  records,  and  frequency 
results  for  the  five  pulse  tests  conducted  are  given  in  Figures  28 
to  52  and  Tables  7  to  11. 


16.2  Smoothing  of  Time  Data 

Inspection  of  the  experimental  time  data  indicates  that  a 
large  amount  of  extraneous  noise  is  superimposed  on  the  temperature 
signal,  the  signal-to-noise  ratio  for  all  five  runs  being  approxi¬ 
mately  four  to  one.  Because  of  the  relatively  constant  amplitude 
and  high  frequency  of  this  extraneous  disturbance  it  was  suspected 
as  being  static  noise,  and  probably  emanating  from  some  electrical 
source  such  as  the  pump  motors  or  lights  in  the  laboratory.  By 
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FIG.  28  TUBE  FLOW  INPUT  PULSE  -  RUN  17 
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FIG.  29  SHELL  OUTLET  TEMPERATURE  RESPONSE  -  RUN  17 
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FIG.  30  NORMALIZED  FREQUENCY  CONTENT  OF  INPUT  PULSE  -  RUN  17 
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FIG.  31  CALCULATED  AMPLITUDE  RATIO  RESPONSE  -  RUN  17 
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FIG.  32  CALCULATED  PHASE  ANGLE  RESPONSE  -  RUN  17 
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FIG.  33  TUBE  FLOW  INPUT  PULSE  -  RUN  18 
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FIG.  34  SHELL  OUTLET  TEMPERATURE  RESPONSE  -  RUN  18 
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FIG.  35  NORMALIZED  FREQUENCY  CONTENT  OF  INPUT  PULSE  -  RUN  18 
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FIG.  36.  CALCULATED  AMPLITUDE  RATIO  RESPONSE  -  RUN  18 
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FIG.  37  CALCULATED  PHASE  ANGLE  RESPONSE  -  RUN  18 
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FIG.  38  TUBE  FLOW  INPUT  PULSE  -  RUN  19 
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FIG.  39  SHELL  OUTLET  TEMPERATURE  RESPONSE  -  RUN  19 


130 


FIG.  40  NORMALIZED  FREQUENCY  CONTENT  OF  INPUT  PULSE  -  RUN  19 
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FIG.  41  CALCULATED  AMPLITUDE  RATIO  RESPONSE  -  RUN  19 
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FIG.  42  CALCULATED  PHASE  ANGLE  RESPONSE  -  RUN  19 
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FIG.  43  TUBE  FLOW  INPUT  PULSE  -  RUN  20 
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FIG.  44  SHELL  OUTLET  TEMPERATURE  RESPONSE  -  RUN  20 
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FIG.  45  NORMALIZED  FREQUENCY  CONTENT  OF  INPUT  PULSE  -  RUN  20 
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FIG.  46  CALCULATED  AMPLITUDE  RATIO  RESPONSE  -  RUN  20 
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FIG.  47  CALCULATED  PHASE  ANGLE  RESPONSE  -  RUN  20 
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FIG.  48  TUBE  FLOW  INPUT  PULSE  -  RUN  21 
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FIG.  49  SHELL  OUTLET  TEMPERATURE  RESPONSE  -  RUN  21 
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FIG.  50  NORMALIZED  FREQUENCY  CONTENT  OF  INPUT  PULSE  -  RUN  21 
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FIG.  51  CALCULATED  AMPLITUDE  RATIO  RESPONSE  -  RUN  21 
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FIG.  52  CALCULATED  PHASE  ANGLE  RESPONSE  -  RUN  21 


TABLE  7 


Experimental  Data  for  Run  17 


Forcing  Function 


Rectangular  pulse  in  tube  fluid 
(9.0  -*  10.7  +  9.0  U.S.  gpm) 


Response  Function 


Flow  Conditions 


Exchanger  Size 


Shell  liquid  outlet  temperature 


Countercurrent,  water-water 


0.75  inches  x  1.5  inches  x  10  feet 


Experimental  Conditions: 


Tube  fluid  flow  (U.S.  gpm)  -  9.0 

Shell  fluid  flow  (U.S.  gpm)  -  20.0 


Tube  fluid  inlet  temperature  (degrees  F)  -  54.1 

Tube  fluid  outlet  temperature  (degrees  F)  -  79.1 

Shell  fluid  inlet  temperature  (degrees  F)  -  155.8 
Shell  fluid  outlet  temperature  (degrees  F)  -  145.5 
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TABLE  8 


Experimental  Data  for  Run  18 


Forcing  Function 


Response  Function 


Flow  Conditions 


Exchanger  Size 


Ramp  pulse  in  tube  flow 
(7.0  -»■  11.0  7.0  U.S.  gpm) 

Shell  liquid  outlet  temperature 

Countercurrent,  water-water 

0.75  inches  x  1.5  inches  x  10  feet 


Experimental  Conditions: 


Tube 

fluid 

flow  (U.S.  gpm) 

7.0 

Shell 

fluid 

flow 

(U.S.  gpm) 

15.0 

Tube 

fluid 

inlet 

temperature 

(degrees 

F) 

53.4 

Tube 

fluid 

outlet 

temperature 

(degrees 

F) 

81.8 

Shell 

fluid 

inlet 

temperature 

(degrees 

F) 

-  163.3 

Shell 

fluid 

outlet  temperature  (degrees  F) 

-  149.9 

'  ■  . 


TABLE  9 


Experimental  Data  for  Run  19 


Forcing  Function 


Triangular  pulse  in  tube  flow 
(8.0  10.5  +  8.0  U.S.  gpm) 


Response  Function 


Flow  Conditions 


Exchanger  Size 


Shell  liquid  outlet  temperature 


Countercurrent,  water-water 


0.75  inches  x  1.5  inches  x  10  feet 


Experimental  Conditions: 


Tube  fluid  flow  (U.S.  gpm)  -  8.0 

Shell  fluid  flow  (U.S.  gpm)  -  20.0 


Tube  fluid  inlet  temperature  (degrees  F)  -  60.1 
Tube  fluid  outlet  temperature  (degrees  F)  -  85.7 


Shell  fluid  inlet  temperature  (degrees  F)  -  157.1 

Shell  fluid  outlet  temperature  (degrees  F)  -  147.2 


X 
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TABLE  10 


Experimental  Data  for  Run  20 


Forcing  Function 


Displaced  cosine  pulse  in  tube  flow 
(8.0  +  10.5  ->  8.0  U.S.  gpm) 


Response  Function 

Flow  Conditions 


Exchanger  Size 


Shell  liquid  outlet  temperature 


Countercurrent,  water-water 


0.75  inches  x  1.5  inches  x  10  feet 


Experimental  Conditions: 


Tube  fluid  flow  (U.S.  gpm)  -  8.0 

Shell  fluid  flow  (U.S.  gpm)  -  15.0 


Tube  fluid  inlet  temperature  (degrees  F)  -  52.9 
Tube  fluid  outlet  temperature  (degrees  F)  -  78.8 


Shell  fluid  inlet  temperature  (degrees  F)  -  160.4 

Shell  fluid  outlet  temperature  (degrees  F)  -  146.5 
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TABLE  11 


Experimental  Data  for  Run  21 


Forcing  Function 


Rectangular  doublet  in  tube  flow 
(8.0  -»■  10.5  8.0  U.S.  gpm) 


Response  Function 


Flow  Conditions 


Exchanger  Size 


Shell  liquid  outlet  temperature 


Countercurrent,  water-water 


0.75  inches  x  1.5  inches  x  10  feet 


Experimental  Conditions: 


Tube  fluid  flow  (U.S.  gpm)  -  8.0 

Shell  fluid  flow  (U.S.  gpm)  -  15.0 

Tube  fluid  inlet  temperature  (degrees  F)  -  52.2 

Tube  fluid  outlet  temperature  (degrees  F)  -  77.5 

Shell  fluid  inlet  temperature  (degrees  F)  -  158.0 

Shell  fluid  outlet  temperature  (degrees  F)  -  144.8 


, 
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simple  observation  of  the  disturbance  it  was  suspected  of  having  a 
zero  mean  average  (z.m.a.).  Since  Fourier  transformation  is  essenti¬ 
ally  an  integration  technique,  the  presence  of  such  high  frequency 
noise  in  the  time  data  would  not  affect  the  amplitude  ratio  and  phase 
angle  results  at  those  frequencies  below  the  noise  level.  Thus,  if  a 
given  set  of  raw  time  data  superimposed  by  high  frequency  noise  with  a 
zero  mean  average  is  Fourier  transformed,  the  frequency  results  should 
be  the  same  as  those  given  by  a  similar  transformation  of  the  time  data 
with  the  noise  smoothed  out.  This  procedure  constitutes  a  simple  test  as 
to  whether  any  noise  observed  in  pulse  testing  time  data  has  any 
deleterious  effects  on  the  final  frequency  results.  The  only  implication, 
however,  is  that  a  smoothing  routine  must  be  used  to  obtain  smooth 
(noiseless)  time  data  from  the  original  record. 

Another  reason  for  smoothing  time  data  is  for  presentation 
purposes.  In  applying  a  smoothing  routine,  however,  the  experimenter 
must  be  careful  not  to  degrade  the  true  response  information  underlying 
the  noise.  For  this  reason  care  should  be  taken  during  smoothing  of 
pulse  testing  data  to  choose  a  technique  that  has  the  ability  to  en¬ 
hance  the  signal- to-noise  ratio  without  marring  the  true  transient 
response.  Although  not  desiring  to  improve  the  appearance  of  the  shell 
outlet  temperature  plots,  smoothing  was  used  in  this  study  to  prove  or 
disprove  the  presence  of  z.m.a.  noise  in  the  temperature  signal. 

Data  may  be  smoothed  digitally  by  many  different  methods  (6) , 
and  no  attempt  is  made  here  to  present  an  in-depth  study  of  such  a 
broad  subject.  The  aim  is  to  outline  a  few  basic  concepts  of  digital 
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smoothing  and  inform  the  reader  of  the  type,  and  the  reasons  for,  the 
smoothing  technique  used  in  this  study. 


One  of  the  simplest  means  of  smoothing  data  is  by  a  moving 
average.  By  this  method,  a  fixed  number  of  consecutive  data  points 
are  considered,  their  ordinates  are  added  together,  and  division  by 
the  number  of  points  is  made  to  obtain  the  average  ordinate  at  the 
centre  abscissa  of  the  group.  This  technique  is  satisfactory  for 
smoothing  data  which  are  assumed  to  be  at  steady  state,  but  tends  to 
degrade  any  time  record  exhibiting  peaks  (for  example,  transient 
response  information  from  a  pulse  testing  experiment).  To  improve 
upon  the  moving  average,  more  emphasis  could  be  placed  upon  the 
central  raw  data  point  than  upon  its  adjacent  points.  Weighting 
each  of  the  raw  data  points  on  the  basis  of  their  proximity  to  the 
central  data  point  is  the  essence  of  a  convolution  operation.  In 
this  procedure,  the  weighting  values  (convoluting  integers)  are 
multiplied  by  their  pertinent  ordinate  value,  and  the  sum  of  the 
products  is  divided  by  a  normalizing  factor  obtained  by  summing  the 
convoluting  integers.  Various  convoluting  functions  are  possible; 
three  symmetrical  ones  are  shown  in  Figure  53  below. 


Moving  Average  Symmetrical  Triangle  Symmetrical  Exponent 

Fig. 53  Various  Symmetr ical  Convoluting  Functions 
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To  use  any  of  the  above  convoluting  functions,  it  is  implicit 
that  the  experimenter  have  time  data  points  available  both  in  the 
past  and  in  the  future  of  the  central  abscissa  values.  In  DDC 
applications,  data  is  usually  smoothed  on-line,  and  any  convoluting 
set,  therefore,  can  only  be  applied  to  data  points  gathered  in  the 
past,  there  being  no  future  data  points  available.  This  type  of 
convoluting  set  is  termed  nonsymmetrical ,  and  is  the  digital  counter 
part  of  hardware  analog  filters.  For  example,  the  nonsymmetrical 
exponential  filter  available  on  the  IBM  1800  has  the  characteristics 
of  an  RC  filter.  For  process  identification  purposes,  this  type  of 
filter  should  be  avoided,  as  it  tends  to  produce  a  first-order  lag 
effect  in  the  smoothed  time  data.  It  is  conceivable  that  compara¬ 
tively  severe  smoothing  of  this  type,  when  applied  to  pulse  testing 
data,  could  cause  significant  errors  in  interpreting  the  resulting 
frequency  information. 

The  convoluting  function  used  in  this  study  is  based  on  the 
concept  of  curve  fitting  by  least  squares.  The  derivation  of  the 
function  is  given  by  Savitzky  and  Golay  (47).  For  either  a  cubic 
or  a  quadratic  function,  the  set  of  convoluting  integers  is  the  same 
and  the  set  for  up  to  25  points  with  the  appropriate  normalizing 
factors  is  given. 

The  experimental  cosine  pulse  data  was  used  to  test  for  the 
presence  of  z.m.a.  noise  overlying  the  temperature  signal.  To 
smooth  the  experimental  temperature  data,  the  11-point  convoluting 
set  given  by  the  above  authors  was  used.  The  raw  and  smooth  data 
were  analyzed;  the  resulting  frequency  information  is  given  in 
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Table  12 .  The  agreement  between  the  results  does  indicate  that  the 


extraneous  disturbance  is,  in  fact,  z.m.a.  noise. 


16.3  Discussion  of  Experimental  Frequency  Results 

Initial  experimentation  indicated  that  a  change  in  the  shell 
fluid  outlet  temperature  resulting  from  a  tube  flow  pulse  would  be 
in  the  order  of  1.5  to  2.0  degrees  F.  Considering  that  normal 
process  fluctuations  result  in  temperature  variations  of  ±  0.25 
degrees  F,  concern  was  felt  that  the  resulting  frequency  information 
would  be  badly  scattered.  On  analyzing  the  data,  however,  the 
frequency  results  were  found  to  be  very  uniform  in  their  appearance 
and  serious  scatter  was  encountered  only  in  the  higher  frequency 
range,  where  it  was  expected.  This  attests  to  the  worth  of  the 
pulse  testing  method. 

The  time  data  sampling  interval  of  0.05  seconds  implies  that 
erroneous  frequency  results  emanating  from  numerical  approximation 
errors  were  to  be  expected  at  frequencies  higher  than  3  radians  per 
second.  Scatter  in  the  frequency  diagrams,  however,  was  apparent  at 
lower  frequencies  than  this  due  to  the  frequency  content  of  the 
pulses  approaching  zero.  It  was  necessary,  therefore,  to  calculate 
the  frequency  content  for  each  pulse  to  determine  the  frequency 
range  of  valid  frequency  information.  As  anticipated,  the  tri¬ 
angular  and  cosine  pulses,  having  shapes  approaching  the  Dirac 
impulse  function,  exhibited  better  frequency  content  characteristics 
than  the  other  one-sided  forcing  functions.  It  is  noticed  in  the 
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TABLE  12 


Comparison  of  Frequency  Results  of 
Raw  and  Smooth  Data  from  Run  20 

Frequency  Results  using  Raw  Data: 


FREQUENCY 
RAD/SEC 
0.0099 
0.0119 
0 . 0  1 A  2 
0.0170 
0.0203 
0.0242 
0.0289 
0  •  0  3  A  5 
0.0A12 
0.049? 
0.0587 
0.0701 
0.0837 
0.0999 
0.1193 
0.1425 
0.1701 
0.2030 
0.2424 
0.2894 
0.3455 
0.4124 
0.4923 
0.5877 
0.7017 
0.8376 
0.9999 
1.1937 
1.4250 
1.7012 
2.0308 
2.4244 
2.8942 
3.4550 
4.1245 
4.9238 
5.8779 
7.0169 
8.3766 
9.9998 


FREQUENCY 

CONTENT 

0.9999 

0.9999 

0.9999 

0.9999 

0.9999 

0.9998 

0.9998 

0.9997 

0.9996 

0.9995 

0.9992 

0.9989 

0.9985 

0.9979 

0.9970 

0.9958 

0.9940 

0.9915 

0.9879 

0.9828 

0.9755 

0.9652 

0.9506 

0.9299 

0.9008 

0.8599 

0.8033 

0.7258 

0.6225 

0.4903 

0.3319 

0.1642 

0.0625 

0.1283 

0.1428 

0.0784 

0.0728 

0.0615 

0.0543 

0.0338 


AMPLITUDE 

RATIO 

1.101 

1.101 

1.100 

1.100 

1.100 

1.099 

1.098 

1.096 

1.094 

1.091 

1.087 

1.081 

1.073 

1.061 

1.044 

1.021 

0.989 

0.946 

0.888 

0.815 

0.725 

0.625 

0.527 

0.450 

0.396 

0.334 

0.226 

0.127 

0.145 

0.001 

0.042 

0.222 

0.631 

0.056 

0.033 

0.311 

0.337 

0.067 

0.470 

0.455 


PHASE 

ANGLE 


-2 

•  362 

-2 

.819 

-3 

.  366 

-4 

.018 

-4 

.796 

-5 

.725 

-6 

.833 

-8 

.156 

-9 

.733 

“11 

.614 

-13 

.857 

“16 

.527 

-19 

.706 

“2  3 

.482 

“27 

.960 

-33 

.253 

-39 

•  481 

-46 

.757 

-55 

.161 

-64 

.694 

-75 

.197 

“86 

.274 

“97 

.505 

-109 

.660 

-126 

.580 

-152 

.314 

-182 

.471 

-184 

.399 

-198 

.880 

83 

.269 

-150 

.782 

-122 

.607 

-181 

.053 

-111 

.CIS 

-267 

.939 

-199 

.  106 

23 

.539 

-101 

.030 

-175 

.779 

39 

.141 

TABLE  12  Continued.... 


Frequency  Results  using  Smooth  Data: 


FREQUENCY 
RAD/SEC 
0.0099 
0.0119 
0  •  0  1 A  2 
0.0170 
0.0203 
0.0242 
0.0289 
0. 0345 
0.0412 
0.0492 
0.0587 
0.0701 
0.0837 
0.0999 
0.1193 
0.1425 
0.1701 
0.203C 
0.2424 
0.2894 
0.3455 
0.4124 
0.4923 
0.5877 
0.7017 
0.8376 
0.9999 
1.1937 
1.4250 
1.7012 
2.0308 
2.4244 
2.8942 
3.4550 
4.1245 
4.9233 
5.8779 
7.0169 
8 .3766 
9.9998 


FREQUENCY 

CONTENT 

0.9999 

0.9999 

0.9999 

0.9999 

0.9999 

0.9998 

0.9998 

0.9997 

0.9996 

0.9995 

0.9992 

0.9939 

0.99S5 

0.9979 

0.9970 

0.9958 

0.9940 

0.9915 

0.9879 

0.9828 

0.9755 

0.9652 

0.9506 

0.9299 

0.9008 

0.6599 

0.8033 

0.7258 

0.6225 

0.4903 

0.3319 

0.1642 

0.0625 

0.1283 

0.1428 

0.0784 

0.0728 

0.0615 

0.0543 

0.0338 


AMPLITUDE 

RATIO 

1.104 

1.104 

1.104 

1.104 

1.103 

1.102 

1.101 

1.100 

1.098 

1.095 

1.091 

1.085 

1.076 

1.064 

1.047 

1.024 

0.992 

0.949 

0.391 

0.817 

0.727 

0.626 

0.528 

0.451 

0.397 

0.333 

0.223 

0.125 

0.145 

0.C04 

0.045 

0.233 

0.612 

0.067 

0.026 

0.323 

0.352 

0.090 

0.449 

0.356 


PHASE 
ANGLE 
-2.371 
-2 .830 
-3.373 
-4.033 
-4.814 
-5.746 
-6.859 
-8.187 
-9.770 
-11.658 
-13.909 
-16.589 
-19.778 
-23.563 
-28.060 
-33.369 
-39.614 
-46.905 
-55.319 
-64.851 
-75.327 
-86.331 
-97.425 
-109.407 
-126.234 
-152.074 
-182.623 
-184.770 
-200. 157 
-253.478 
-155.373 
-122.691 
-182.243 
-119.468 
77.488 
-20 1 . 644 
23.838 
-102.831 
-172.644 
35.556 


frequency  results  that  serious  scatter  and  obviously  erroneous 
amplitude  ratio  and  phase  angle  information  is  exhibited  at  those 
frequencies  for  which  the  corresponding  normalized  frequency  content 
curve  is  less  than  about  0.3. 

The  rectangular  doublet  pulse  gives  very  poor  steady  state 
gain  information.  This  is  because  the  amplitudes  of  the  component 
frequencies  of  this  type  of  pulse  approach  zero  at  the  lower 
frequencies.  From  its  frequency  content  curve,  the  component 
frequency  of  the  pulse  exhibiting  a  maximum  amplitude  occurs  at 
approximately  1.5  radians  per  second.  This  type  of  pulse  which 
exhibits  a  maximum  in  its  frequency  content  curve,  although  reveal¬ 
ing  very  poor  steady  state  gain  information,  could  be  very  useful 
for  studying  desired  frequency  ranges.  For  example,  the  phenomenon 
of  "resonance",  explained  in  detail  by  Fisher  (23)>  occurs  in  heat 
exchangers  at  a  frequency  of  2ttV/L  radians  per  second,  where  L  is 
the  length  of  the  heat  exchanger  in  feet,  and  V  is  the  velocity  in 
feet  per  second,  of  the  fluid  whose  temperature  is  being  measured. 
For  the  conditions  of  this  study,  the  resonant  frequency  of  the  test 
heat  exchanger  would  occur  at  approximately  3  radians  per  second. 
This  is  unfortunately  at  the  limit  of  the  valid  frequency  range  of 
the  rectangular  doublet  pulse  run,  and  thus  the  resonance  phenomenon 
is  not  apparent  in  the  resulting  frequency  information.  Although 
an  investigation  of  resonance  is  beyond  this  study,  any  future  work 
should  involve  shifting  the  range  of  valid  frequency  information  to 
encompass  the  frequency  of  3  radians  per  second.  This  would  entail 
using  all  40  points  per  second  in  the  time  data  to  ensure  that  any 
numerical  approximation  problems  would  be  obviated,  and  shortening 


■ 
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the  width  of  the  rectangular  doubletpulse  from  4  to  2  seconds  to  shift 
the  maximum  frequency  content  to  a  higher  frequency. 

In  the  amplitude  ratio  and  phase  angle  plots  for  all  the  pulse 
tests,  a  small  peak  is  seen  to  occur  at  a  frequency  ranging  between  0.5 
and  1.0  radians  per  second.  The  peak  occurs  if  the  experimental  data  is 
analyzed  using  a  random  noise  analysis  program,  indicating  its  existence 
is  due  to  some  physical  phenomenon  rather  than  the  method  of  analysis. 

This  conclusion  is  also  supported  by  the  fact  that  it  appears  in  the 
data  of  Lees  (35),  the  author's  predecessor  in  this  work.  Although  its 
occurrence  is  not  understood,  it  is  possible  that  the  peak  could  be  due 
to  the  change  in  turbulence  conditions,  and  thus  heat  transfer  rates, 
which  result  during  the  implementation  of  the  pulse.  Another  possibility 
is  that  it  may  be  due  to  the  dynamics  of  the  tube  and  shell  walls .  This 
is  a  matter  of  conjecture,  however,  and  should  be  investigated  further. 


16.4  Simple  Heat  Exchanger  Models 

The  general  four-equation  model  describing  the  dynamics  of  a  heat 
exchanger,  and  used  in  Part  I  of  this  thesis,  would  probably  not  find 
industrial  acceptance  because  of  its  complexity.  For  engineering 
design  purposes,  much  simpler  models  are  desired. 

Work  done  by  Morris  (42) ,  Hearn  (26)  and  Vincent  (56)  indi¬ 
cated  that  a  complex  nonlinear  heat  exchanger  may  be  approximated 
with  reasonable  success  by  a  simple  first-order  equation.  Inspection 
of  the  transient  responses  obtained  in  this  study  suggests  that  the 
approximation  would  be  improved  by  using  a  transportation  lag  term 
in  series  with  the  first-order  system.  The  transfer  function  of  such 
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a  model  (Model  1)  may  therefore  be  expressed  as: 

G(s)  =  e~T2S  ...  16-1 

(1  +  tj  s) 

where  G(s)  is  the  transfer  function  relating  the  shell  fluid  outlet 
temperature  to  a  tube  flow  disturbance. 

From  the  frequency  results  of  a  pulse  test,  it  is  possible 
to  estimate  values  for  the  time  constants  x1  and  for  Model  1. 

Consider  the  characteristics  of  the  two  contributing  systems.  For 
the  first-order  system: 


A.R.  =  1  ...  16-2 

/x2  oo2  +  1 

1 


<J>  =  tan  (-  toil)  degrees  ...  16-3 

and  for  the  transportation  lag: 


A. R.  =1.0  ...  16-4 

<p  =  -  L0T2  radians  •  •  •  16-5 

The  transportation  lag,  having  an  amplitude  ratio  of  unity, 
means  that  the  amplitude  ratio  characteristics  of  Model  1  are  the 
same  as  those  for  the  first-order  system  alone.  The  time  constant, 
x 1 ,  is  calculated  from  the  reciprocal  of  the  corner  frequency,  the 
corner  frequency  occurring,  of  course,  at  an  amplitude  ratio  equal 
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to  0.707  of  the  steady  state  gain.  At  the  corner  frequency,  the 
phase  lag  contributed  by  the  first-order  dynamics  of  the  process 
equals  45  degrees.  The  phase  lag  value  over  and  above  this  is 
therefore  considered  to  be  due  to  the  transportation  lag.  Dividing 
this  phase  lag  difference,  expressed  in  radians,  by  the  corner 
frequency,  gives  a  value  for  t2. 


Calculation  of  the  time  constants  for  Model  1  is  based  on 
the  premise  that  an  experimental  Bode  diagram  is  available.  To 
circumvent  this  problem,  Thal-Larsen  (55)  presents  a  model  in  which 
the  time  constants  are  obtained  from  knowledge  of  the  steady  state 
operation  of  a  heat  exchanger.  In  developing  the  model,  the 
experimental  frequency  response  data  of  a  number  of  heat  exchangers, 
including  both  liquid-liquid  exchangers  and  condensers,  obtained  by 
other  workers  were  utilized.  Thal-Larsen  showed  that  the  normalized 
frequency  response  diagrams  were  very  similar,  and  further  demon¬ 
strated  that  the  experimental  frequency  plots  could  be  approximated 
by  a  simple  model  for  which  the  time  constants  are  only  functions  of 
fluid  residence  times.  The  model,  for  the  purposes  of  this  thesis, 
will  be  called  Model  2. 


G(s)  = 


-ToS 


(TjS  +  1)(t2S  +  l) 


,  R  +  R 

where  ii  =  A  1 


Ra  +  Ri 

T  2  =  A  1 


. ..  16-6 


8 


' 
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and  R.  =  fluid  residence  time  in  the  annulus 
A 

=  fluid  residence  time  in  the  tube. 

The  time  constants  for  both  models  were  calculated  for  the  one 
sided  pulse  tests  conducted  in  this  study;  the  values  appear  in 
Table  13.  The  experimental  frequency  response  and  those  predicted 
by  Models  1  and  2  are  shown  in  Figures  54  to  63. 

Considering  the  simplicity  of  the  models,  both  give  very 
reasonable  approximations  to  the  experimental  curves.  Particularly 
edifying  is  the  agreement  between  Model  2  and  the  experimental 
results  as  it  indicates  a  very  convenient  engineering  approach  for 
predicting  the  dynamics  of  heat  exchangers. 


Time  Constants  for  Model  1  and  Model 
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FIG.  54  COMPARISON  OF  THE  AMPLITUDE  RATIO  PREDICTIONS 
AS  GIVEN  BY  MODEL  1  AND  MODEL  2  FOR  RUN  17 
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FIG.  55  COMPARISON  OF  THE  PHASE  ANGLE  PREDICTIONS 
AS  GIVEN  BY  MODEL  1  AND  MODEL  2  FOR  RUN  17 
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FIG.  56  COMPARISON  OF  THE  AMPLITUDE  RATIO  PREDICTIONS 
AS  GIVEN  BY  MODEL  1  AND  MODEL  2  FOR  RUN  18 
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FIG.  57  COMPARISON  OF  THE  PHASE  ANGLE  PREDICTIONS 
AS  GIVEN  BY  MODEL  1  AND  MODEL  2  FOR  RUN  18 
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FIG.  58  COMPARISON  OF  THE  AMPLITUDE  RATIO  PREDICTIONS 
AS  GIVEN  BY  MODEL  1  AND  MODEL  2  FOR  RUN  19 
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FIG.  59  COMPARISON  OF  THE  PHASE  ANGLE  PREDICTIONS 
AS  GIVEN  BY  MODEL  1  AND  MODEL  2  FOR  RUN  19 
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FIG.  60  COMPARISON  OF  THE  AMPLITUDE  RATIO  PREDICTIONS 
AS  GIVEN  BY  MODEL  1  AND  MODEL  2  FOR  RUN  20 
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FIG.  61  COMPARISON  OF  THE  PHASE  ANGLE  PREDICTIONS 
AS  GIVEN  BY  MODEL  1  AND  MODEL  2  FOR  RUN  20 
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17  CONCLUSIONS 

A  successful  user-orientated  computerized  system  for  on-line  pulse 
testing  was  designed  and  demonstrated. 

A  simple  straightforward  set  of  guidelines  (Chapter  13)  for  designing 
a  pulse  testing  experiment  has  been  presented  in  the  hope  that  it 
might  promote  more  widespread  use  and  acceptance  of  this  process 
identification  technique  in  industry. 

For  any  industrial  application,  probably  the  main  concern  of  any 
pulse  testing  experiment  is  that  it  should  reveal  satisfactory  steady 
state  gain  and  low  frequency  information.  This  would  suggest  that 
one-sided  pulses  are  to  be  preferred  in  plant  testing.  Only  in 
special  cases,  for  example  in  resonance  studies  on  heat  exchangers 
where  a  band  of  higher  frequencies  is  of  interest,  should  double¬ 
sided  pulses  be  used. 

From  frequency  content  considerations,  sharp-pointed  symmetrical 
pulses  are  favoured  for  pulse  testing.  The  experimental  frequency 
results  obtained  from  the  one-sided  pulse  runs  conducted  in  this 
study,  however,  show  that  the  difference  in  regard  to  the  apparent 
range  of  valid  frequency  data  is  quite  marginal,  scatter  in  the  fre¬ 
quency  results  beginning  in  all  cases  at  a  frequency  of  approximately 
1.2  radians  per  second.  This  would  suggest  that  in  any  industrial 
testing  no  concern  should  be  made  about  attempting  to  realize  the 
"better"  pulse  shapes. 


5. 


The  experimental  frequency  results  appear  to  be  valid  down  to  a 
normalized  frequency  content  of  approximately  0.4;  this  compares 
with  the  figure  of  0.3  which  is  usually  cited  in  the  literature. 
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This  demonstrates  that  no  specific  rule  can  be  made  about  the 
relation  of  valid  frequency  data  and  normalized  frequency  content. 

As  already  indicated  in  Chapter  10,  there  are  numerous  factors 
affecting  the  accuracy  of  pulse  testing  frequency  results,  and 
although  the  frequency  content  of  the  input  pulse  may  be  a  major 
factor  it  is  not  the  only  one.  Fortunately  for  the  experimenter, 
invalid  frequency  data  are  easily  recognizable  as  "scatter"  in  the 
frequency  diagrams. 

6.  In  Chapter  11  a  comparison  of  data  analysis  methods  was  made,  and 

one  aspect  discussed  was  the  relative  computing  speeds  of  the  methods. 

To  ensure  that  this  factor  is  kept  in  context  it  should  be  mentioned 
that  the  computer  costs  associated  with  a  pulse  testing  experiment 

Probably  be  negligible  in  comparison  to  the  associated  plant  costs. 

7.  Any  time  data  intended  for  dynamic  analysis  studies  should,  if  possible, 
be  collected  with  the  minimum  of  smoothing  (digital  and  analog)  in  order 
to  avoid  artifically  induced  time  lags  being  exhibited  in  the  smoothed 
data. 

8.  Extraneous  noise  overlying  transient  response  information,  if  from  an 
electrical  source,  will  probably  not  affect  any  Fourier  analysis 
results,  at  least  in  the  frequency  region  of  concern.  If,  as  a  test 
for  zero  mean  average  noise  or  for  presentation  purposes,  it  is 
desired  to  smooth  any  pulse  testing  time  data,  then  a  symmetrical 
convoluting  function  of  the  type  described  in  this  study  should  be  used. 

9.  If  a  simple  model  of  a  heat  exchanger  is  desired  for  control  design 
purposes  or  as  an  aid  in  designing  a  pulse  testing  experiment,  the 
model  given  by  Thal-Larsen  and  designated  Model  2  in  this  study 


should  be  satisfactory. 
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18  RECOMMENDATIONS  FOR  FUTURE  WORK 


Alteration  to  Equipment 

In  the  experimental  study,  the  magnitudes  of  the  temperature 
responses  resulting  from  the  tube  fluid  pulses  were  of  the 
order  of  1.0  to  1.5  degrees  F.  Temperature  changes  due  to 
natural  process  fluctuations  were  of  the  order  of  0.5  degrees 
F,  thus  making  it  sometimes  difficult  to  recover  useful  time 
data.  For  future  pulse  testing  studies  it  is  recommended 
that  the  heat  exchanger  tube  (I.D.  =  0.75  inches)  be  replaced 
by  one  having  a  larger  internal  diameter,  say,  one  inch.  The 
resulting  change  in  the  ratio  of  the  volumetric  flow  rates 
occurring  in  the  annulus  and  tube  should  produce  a  higher 
gain  system;  i.e.,  a  higher  shell  fluid  outlet  temperature 
for  a  given  change  in  the  controller  output  to  the  tube  fluid 
control  valve.  This  would  tend  to  minimize  errors  in  the 
frequency  diagrams  caused  by  any  uncontrolled  process  dis¬ 
turbances  . 

Improvement  of  the  High  Speed  Data  Acquisition  Program  (VCDAQ) 

The  present  high  speed  data  acquisition  program  is  written 
in  such  a  manner  that,  at  the  commencement  of  a  pulse 
testing  run,  the  DDC  function  of  the  computer  is  suspended. 

Any  extraneous  disturbances  occurring  in  the  process  during 
this  period  tend,  therefore,  to  invalidate  the  pulse  testing 
data.  The  program  VCDAQ  should  be  modified  to  correct  this 


problem. 
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Heat  Exchanger  Transfer  Functions 

This  study  was  only  concerned  with  the  transfer  function 
relating  shell  fluid  outlet  temperature  to  tube  fluid  flow. 
Consideration  should  be  given  to  program  and  equipment 
changes  necessary  to  obtain  other  transfer  functions  of 
the  heat  exchanger. 


4.  CSMP  Model 

Using  the  CSMP  model  developed  in  Part  I  of  this  thesis  the 
experimental  pulse  testing  runs  should  be  simulated.  The 
resulting  time  data  should  be  analyzed  and  the  frequency 
results  compared  with  those  given  by  experiment. 

5 .  The  Peak  in  the  Frequency  Plots 

An  investigation  should  be  made  into  the  cause  of  the  small 
peak  observed  in  the  Bode  diagrams,  occurring  at  a  frequency 
of  between  0.5  and  1.0. 

6.  Resonance 

Advantage  should  be  taken  of  the  peculiar  frequency  content 
curve  of  the  double-sided  pulse  (i.e.,  that  it  exhibits  a 
maximum  at  some  frequency  greater  than  zero) ,  in  a  study  of 


resonance. 
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Modeling  Studies 


Further  experimental  studies  should  yield  amplitude  ratio  and 
phase  angle  information  which  is  accurate  to  higher  frequencies 
than  obtained  in  this  work.  Availability  of  such  information 
would  allow  the  study  of  more  rigorous  models  proposed  in  the 
literature  to  describe  heat  exchanger  dynamics.  Of  particular 
interest  would  be  the  model  given  by  Fisher  (23),  due  to  its 
ability  to  predict  resonance. 

Random  Noise  Analysis 

Pulse  testing  requires  the  use  of  an  input  forcing  function 
which  has  to  be  generated  and  purposely  implemented  by  the 
user.  An  indirect  method  of  process  characterization,  random 
noise  analysis,  utilizes  the  records  of  the  inputs  and  outputs 
obtained  in  the  normal  operation  of  the  system  (36) .  This 
has  obvious  advantages.  A  comparative  study  of  the  random 
noise  and  pulse  testing  analysis  methods  is  recommended. 
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NOMENCLATURE 

General 

A 

cross  sectional  area  -  ft2 

Cp 

heat  capacity  -  Btu/lb  °F 

G(  )  = 

transfer  function 

h 

heat  transfer  film  coefficient  -  Btu/ft2  sec  °F 

• 

J 

square  root  of  -1 

k 

thermal  conductivity  -  Btu/sec  ft  °F 

L 

heat  exchanger  length  -  ft 

P 

perimeter  -  ft 

Q 

heat  transferred  -  Btu/sec 

s  = 

Laplace  transform  operator  -  sec  ^ 

t 

time  -  sec 

T 

temperature  -  °F 

U 

overall  heat  transfer  film  coefficient  -  Btu/sec  ft2  °F 

v  = 

velocity  -  ft/sec 

x  = 

length  -  ft 

O)  = 

frequency  -  radians/sec 

p 

density  -  lb/ft3 

y 

viscosity  -  lb/ft  sec 

i  = 

time  constant  -  sec 

Subscripts 


1 

tube  fluid 

2 

annular  fluid 

S 

shell 

T 

tube 

.  -  o  3ocn  J  • 
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APPENDIX  I-A 


Special  CSMP  Subroutines 


In  the  following  listing  there  are  four  CSMP  subroutines  which 
constitute  special  elements  .  The  function  of  these  elements  is 
as  follows: 

(i)  Subroutine  2  generates  a  variable  period,  variable  magnitude 
rectangular  pulse: 


where  Parameter  1  =  time  before  initialization  of  the  pulse 
Parameter  2  =  period  of  the  pulse 
Parameter  3  =  magnitude  of  the  pulse 

Note  that,  for  this  element,  during  the  nonpulsed  period,  the 
output  is  equivalent  to  the  input. 

(ii)  Subroutine  3  generates  a  tube  film  coefficient  as  a  function 
of  liquid  bulk  temperature  and  velocity.  The  relationship 
is  given  by  the  Dittus-Boelter  equation, 

4id\  =  0.023  (Re)°‘8  (Pr)°*4 
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where  the  subscript  b  denotes  that  the  properties  are  to  be 
evaluated  at  the  average  bulk  temperature  of  the  stream. 


Parameter  1 
Input  1 
Input  2 
Output 


Tube  I.D.  (inches) 

Tube  fluid  velocity  (ft/sec) 
Temperature  (degrees  F) 
h  (Btu/sec  ft2  degrees  F) 


(iii)  Subroutine  4  generates  a  tube-to-annulus  film  coefficient  based 
on  McAdams'  modification  of  the  Dittus-Boelter  equation: 


In  the  program, 

Parameter  1 
Parameter  2 
Input  1 
Input  2 
Output 


Shell  I.D.  (inches) 

Tube  O.D.  (inches) 

Average  fluid  velocity  (ft/sec) 
Average  fluid  temperature  (degrees  F) 
h  (Btu/sec  ft2  degrees  F) 


(iv)  Subroutine  5  generates  a  shell  film  coefficient  using  the 
equation 
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In  this  subroutine, 


Parameter  1  = 

Shell  I.D.  (inches) 

Parameter  2  E 

Tube  O.D.  (inches) 

Input  1  E 

Velocity  (ft/sec) 

Input  2  e 

Average  bulk  temperature  of  annular  fluid 
(degrees  F) 

Output  E 

h  (Btu/sec  ft2  degrees  F) 

(  £  V:  a  i! 
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SUBROUTINE  SUB2 
SUBROUTINE  SUB2 


this  element  generates  a  variable  period*variable 

MAGNITUDE  PULSE 

PAR  1 (  I ) c  T I  ME  BEFORE  INITIALIZATION  OF  THE  PULSE 
PAR2 ( I ) *PER I OD  OF  PULSE 
PAR3(I)=  MAGNITUDE  OF  PULSE 


REAL  RE ALS ( 395 ) 

INTEGER  I  NTS ( 58 7 ) 

DIMENSION  C(76)*MTRX2(75) *MTRX3 ( 75 ) *MTRX4 ( 75 ) 
DIMENSION  PARI ( 75 ) *PAR2 ( 75 ) *PAR3 (75) 

COMMON  REALS  *  I  NTS 
EQUIVALENCE  (REALS  (2) 

( REALS (81) ♦ 

(REALS! 156 ) 

( REALS ( 231 ) 

( REALS ( 77 ) 


EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 


(  I  NTS ( 76 ) 

( INTS( 151 ) 
( I  NTS ( 226 ) 
( I  NTS ( 376 ) 


♦  C  ( 1 )  ) 

PARI ( 1 ) ) 
PAR2 ( 1 ) ) 
PAR3 ( 1 ) ) 

T  ) 

MTRX2 ( 1 ) ) 
MTRX3 ( 1 ) ) 
MTRX4 ( 1 ) ) 
I  ) 


I  CONTAINS  THE  INDEX  FOR  THE  OUTPUT  VARIABLE 

MTRX2 (  I  )  CONTAINS  THE  BLOCK  NUMBER  *  J  »  OF  THE  l‘ST 

INPUT  TO  BLOC 

C(J)  CONTAINS  THE  CURRENT  VALUE  OF  THE  1»ST  INPUT  TO 
BLOCK  I 

THE  CURRENT  VALUE  OF  THE  OUTPUT  FROM  BLOCK  I  IS  STORED 
IN  C ( I  ) 


J=MTRX2 ( I ) 
C J=C ( J ) 


C  CALCULATIONS 


IF  ( PAR  1 ( I ) -T  >  2*2*1 

1  C(I)=CJ 
GO  TO  5 

2  IF  (PARI (  I  J+PAR2 ( I )-T )  4*3*3 

3  C ( I ) =PAR3 ( I ) *CJ 
GO  TO  5 

4  C(I)=CJ 

5  RETURN 
END 
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SUBROUTINE  SUB3 
SUBROUTINE  SUB3 

THIS  SUBROUTINE  CALCULATES  THE  INSIDE  TUBE  FILM 
COEFFICIENT 


IN  THIS  CALC  H  =  F ( VELOC I TY *  AVERAGE  BULK  TEMP) 
LET  C( I )aH=BTU/SEC»FT«**2«Q»DEG  F 
LET  C ( J)=VELOCITY  ( FT/SEC ) =1 • ST  INPUT 
LET  C ( K ) ^TEMPERATURE  ( DEG«  F ) =2 1 ND  INPUT 
LET  PAR  1 (  I )  =  TUBE  I*D.  (INCHES) 


REAL  REALS ( 395 ) 

INTEGER  I  NTS ( 587 ) 

DIMENSION  C ( 76 ) *MTRX2(75) *MTRX3(75) *MTRX4(75) 
DIMENSION  PARI ( 75 ) *PAR2 (75) *PAR3 ( 75 ) 

COMMON  REALS* INTS 
EQUIVALENCE  (REALS  (2)  *  C ( 1 ) ) 

EQUIVALENCE  ( REALS ( 81 ) *  PARI ( 1 ) ) 

EQUIVALENCE  ( REALS ( 156 ) *PAR2 ( 1 ) ) 

EQUIVALENCE  ( REALS ( 231 ) *  PAR3 ( 1 ) ) 

EQUIVALENCE  (INTS(76)  *  MTRX2 ( 1 ) ) 

EQUIVALENCE  < INTS ( 151 )  #  MTRX3 ( 1 ) ) 

EQUIVALENCE  ( I  NTS  <  226 )  *  MTRX4 ( 1 ) ) 

EQUIVALENCE  (INTS(376)  *  I  ) 

I  CONTAINS  THE  INDEX  FOR  THE  OUTPUT  VARIABLE 
C(J)  CONTAINS  THE  CURRENT  VALUE  OF  THE  l’ST  INPUT  TO 
BLOCK  I 

C ( K )  CONTAINS  THE  CURRENT  VALUE  OF  THE  2 • ND  INPUT  TO 
BLOCK  I 

J=MTRX2 ( I ) 

K=MTRX3 ( I ) 

CALCULAT IONS 

CONVERT  TEMP  FOR  VISCOSITY  CALCULATION 

TC=(C(K>-32#0)*(5.0/9#0> 

C  CALCULATE  VISCOSITY 

VIS=( TC-8. 435) +SQRT( 8078. 4+( ( TC-8 .435 ) **2 . 0 ) ) 

VIS=2»1482*VIS-120«Q 

VIS=( 1.0/VIS)*0.0672 


C 


CALCULATE  THE  RE • NO« 
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SUBROUTINE  SUB3 

RE= ( C ( J ) *PAR 1 (I)*62«4)/(VIS*12«0) 

C  CALCULATE  THERMAL  CONDUCTIVITY 

1*C?KH (0#00266*( C (K ,**3*0 > )-<3.2*<C(K>**2.0> )+( 1073.33 
8+286000*0 ) / 10000 00*0 

C  CALCULATE  THE  PRANDTL  NO. 

PR= ( VIS*3600.0) /TCON 

C  CALCULATE  THE  NUSSELT  NO. 

ANU=0 .023* ( RE**0 •8)*(PR**0.4) 

C  CALCULATE  THE  FILM  COEFFICIENT 

C ( I )  =  ( ANU*TC0N*12 .0 ) / ( PARI ( I) *36  00.0  ) 

RETURN 

END 
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SUBROUTINE  SUB4 
SUBROUTINE  SUB4 


THIS  SUBROUTINE  CALCULATES  THE  OUTSIDE  TUBE  FILM 
COEFFICIENT 


H=F< VELOCITY i AVG.BULK  TEMPERATURE) 

LET  C ( I )=H=BTU/SEC.FT***2.0.DEG  F 
LET  C ( J ) = VELOC I T Y  ( FT • /SEC )  =  1  *  ST  INPUT 
LET  C  (  K  )  =  AVG  •  FLUID  TEMP  (  DEG*F  ) *2 • ND  INPUT 
LETC(L)=AVG.  WALL  TEMP.  ( DEG* F ) =3 1 RD  INPUT 
LET  PAR  1  (  I  ) =SHELL  I.D*  (INCHES) 

LET  PAR2 (  I  ) -TUBE  O.D.  (INCHES) 

LET  PAR 3 (  I  )=CORRELATIONCOEFFICIENT 


REAL  RE ALS ( 395 ) 

INTEGER  I  NTS ( 587 ) 

DIMENSION  C ( 76 ) *MTRX2 ( 75 ) #MTRX3 ( 75 ) *MTRX4 ( 75 ) 
DIMENSION  PARI (75) *PAR2 < 75 ) *PAR3 ( 75 ) 

COMMON  REALS* INTS 
EQUIVALENCE  (REALS  (2)  *  C (1)) 

( REALS (81) *  PARI ( 1 ) ) 

(REALS( 156 ) »  PAR2 ( 1 ) ) 

( REALS ( 231 ) •  PAR3 ( 1 ) ) 

( I  NTS ( 76 ) 


EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

J=MTRX2 ( I ) 

K=MTRX3 ( I ) 

L=MTRX4 ( I ) 


( INTS ( 151 ) 
( INTS ( 226 ) 
( I  NTS ( 376 ) 


*  MTRX2 ( 1 ) ) 

*  MTRX3 ( 1 ) ) 

*  MTRX4 ( 1 ) ) 
»  I  ) 


C  CALCULATIONS 


S I D=PAR 1 (  I  ) 

TOD  =  PAR2 (  I  ) 

C  LET  TF  EQUAL  THE  BULK  SHELL  FLUID  TEMPERATURE 

TF=C ( K ) 

C  CONVERT  TEMP  FOR  VISCOSITY  CALCULATION 

TC=(TF-32.0)*(5. 0/9.0) 

C  CALCULATE  VISCOSITY 

VIS= ( TC-8.435 )+SQRT ( 8078. 4+ ( ( TC-8.435 ) **2.0 ) ) 
VIS=2.1482*VIS-120.Q 
V I S  =  < 1.0/VIS) *0*0672 
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SUBROUTINE  SUB4 

CALCULATE  THE  EQUIVALENT  DIAMETER 

EQD=<SID-T0D)/12.0 

CALCULATE  THE  RE • NO* 

RE= ( C ( J ) *EQD*62.4 ) /VIS 


1  •  A 


CALCULATE  THERMAL  CONDUCTIVITY 

TCON= (( O* 00266* ( TF  **3 • 0  )  ) - ( 3 • 2* ( TF  **2 • 0 ) ) + ( 10 73 • 3 3 
i*tf  )  * 

B+286000 • 0 ) / 10000 00.0 
CALCULATE  THE  PRANDTL  NO. 

PR=(VIS*3600.0) / T  CON 
CALCULATE  THE  NUSSELT  NO. 

ANU=Q.023*(RE**0.8)*(PR**Q.3)*( ( S I D/ TOD >  **0 .45 ) 

CALCULATE  THE  FILM  COEFFICIENT 

C(  I  )  =  (ANU*TCON) / <EQD*3600.0) 

RETURN 

END 


w iTuoaaue 


flJTiviAlCj  7  W  3  J  A  V  i  'JO  3  3f-T  3  T  A  JIQJAO 

o*s:\ (aoT-aia ) -gob 

*OM.3fc  3HT  3T A JUDJA2 
c  1  v\  <  A  *  Sid-*  JOB  *  IU)  )  -3fl 
YT IVI T3UGK0D  JAN'flSHT  3T ABUDJA J 
£  £  ,  )  X  )  +  (  l  0  •  ->  *  *  ■'■  7  )  '••-  S  •  c.  )  •  l  (  G  •  < 

(  "i  T  ■'■■  X 
,  .  jL  \  (  C.OOOdSS+B 

♦  Olrt  jra^A^M  BHT  BTAJ'JB.JA^ 

DT \ <O.OOd£*dIV>«^ 

.V,  I  Jjeau/  3hT  3T ABUDJA} 

^  \ 

( Q O T  \Q 1 6  1  )  *  (  £  •  0  *  *  ^  )  *i8»U**3fl)*£SU«C- 

T1/131DI  33  JOD  MJ I 3  3HT  3TAJUDJAO 

(  0  •  G 0  ' fc‘  1  ..  4 )  \  i  ;  ’OT  * Ui'.  A  5  -  ■  I  )  D 

/HUT3F. 

CM3 


nr>  onr»r»o  r»nnnnonn 


1  •  A 


6 


SUBROUTINE  SUB5 
SUBROUTINE  SUB5 

THIS  SUBROUTINE  CALCULATES  THE  INSIDE  SHELL  FILM 
COEFFICIENT 

LET  C (I ) =H=BTU/SEC.FT.**2  .O.DEG.F 
LET  C ( J)=VELOCITY  ( FT/SEC ) =1 » ST  INPUT 
LET  C ( K )= TEMPERA TURE  ( DEG • F • ) =2 ' ND  INPUT 
LET  PARI (  I  )  =SHELL  I.D.  (INCHES) 

LET  PAR2 (  I  )  =  TUBE  O.D.  (INCHES) 


REAL  RE ALS ( 395 ) 

INTEGER  I  NTS ( 58 7  ) 

DIMENSION  C ( 76 ) *MTRX2 ( 75 ) *MTRX3 ( 75 ) *MTRX4 ( 75 ) 
DIMENSION  PARI ( 75 ) *PAR2 (75)*PAR3(75) 

COMMON  REALS* INTS 
EQUIVALENCE  (REALS  (2)  *  C(l)> 

EQUIVALENCE  ( REALS ( 81 ) *  PARI ( 1 ) ) 

EQUIVALENCE  ( REALS ( 1 56 )# PAR2 ( 1 ) ) 

EQUIVALENCE  ( REALS ( 231 ) *  PAR3 ( 1 ) ) 

EQUIVALENCE  (INTS(76)  ♦  MTRX2 ( 1 ) ) 

EQUIVALENCE  (INTS(151)  .  MTRX3 ( 1 ) ) 

EQUIVALENCE  (INTS(226)  ♦  MTRX4 ( 1 ) ) 

EQUIVALENCE  (INTS(376)  *  I  ) 

I  CONTAINS  THE  INDEX  FOR  THE  OUTPUT  VARIABLE 


C(J)  CONTAINS 
BLOCK  I 

THE 

CURRENT 

VALUE 

OF 

THE 

1  ‘ST 

INPUT 

TO 

C ( K )  CONTAINS 
BLOCK  I 

THE 

CURRENT 

VALUE 

OF 

THE 

2  1  ND 

INPUT 

TO 

J=MTRX2 ( I ) 
K=MTRX3 ( I ) 


CALCULATIONS 
EQUIVALENT  DIAMETER 

EQD=PAR 1 ( I J-PAR2 ( I ) 

C  CONVERT  TEMP  FOR  VISCOSITY  CALCULATION 

TC=(C(K)-32#0)*(5. 0/9.0) 

C  CALCULATE  VISCOSITY 

VIS=( TC-8.435 )+SQRT(8078.4+( ( TC-8.435 ) **2.0 ) ) 
VI Ss2.1482*VIS- 120.0 
VIS=( 1 .O/VIS) *0.0672 


C 


CALCULATE  THE  RE. NO. 


edua  ;-wi  rucpaua 


A  •  . 


58ue  3/iTuoflaua 

v.ji3  jjd-,a  icna/.i  3ht  83TA JU3JA3  3imuofl8ua  eiHT 

T/ M3I^H3Q3 

...  .  .  T 

;  TS,X*03fc\T^>  Yin; -i'JV-  U-  >  '  3  i 

QK»  £*  (  *3 #030)  IflU  3T*t*)0  T3J 

(£iH>lI!  •a#  I  JJ3H8*  (  I  )  I  T3J 
(83H3HI)  .r.O  3eui  =  tl)SWAq  7  3J 


te?£)eJA3«  JA3fl 
(V8e>£THl  S3C3TtfI 
[r  T‘  *  {  £Y  )  EXHTM*  (  5V  )  iiX  TM*  ( <;Y  )3  /QI8MMIC 

(e*;*>6«Aq.  iev  )s  *q«(eui*Aq  noianBMia 

sT/i#aJA3fl  kommo3 
((1)3  *  (S)  8JA3H)  'D/13JAVIUG3 

{  (  X  )I>!Aq  •  (  1 8  )  e  J  A 3 ^  JAVIU03 

(  [  i  )S:^Aq.  (  i  51  )d  JA3fl>  33/3  JAV  I  UQ3 
(U)€MAq  *(irSvaJA'^q)  33tf3JAVIU03 
IDSXHTM  ,  (  d  V  )  £  T  kl  I  )  23iA3JAVIUQ3 

taitxRTM  ♦  (xeneTkin  3Di*3javiuo3 
((i)Ax;;iTi  *  '  ?.<s)aTtfi)  o/3jAviud3 

(  I  «  (dU)2Tk',!>  :OH3  J  AV I UOl 
» 

J  f  -  A  I  fl  AV  TU9TUOr3HT  >*03  X3Cini  3HT  81*11  AT /CO  I 
.  .  i  ,  ,  .  .  •  O  3k  '  i  A  .  DD  (U^ 

i  >30ja 

;  3HT  eMJATWOD  {>>3 

I  X30JQ 

{  I  )  SX^TY.rc 
{  I  )  EXflTM=X 

8/01  TAJU3JA3 
«3T3MAia  T/.3JAVIUQ3 

( i )  $v*q-(  i )  i«Aq=*ao3 

HO  I TA JLO  J A3  YTI80D8IV  WC3  3M3T  T83V/03 

{ (- . e\0 •  )J*  (  .££-(>)  3  JOT 

YT I £038  I  V  3TAJU3JA0 


.  .  .  .  -  .  .  (eC^«8-3T) » 2 1 V 

C.  i  ~  c  1  V  "•■ s  3AX*S  =  eiV 

sva. .o* i a i v \ o . i  ) =e i v 

3HT  3TAJU3JA0 


1  •  A 
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SUBROUTINE  SUB5 

RE-(C(J)*EQD  *62.4)/ (VIS*12.0) 

C  CALCULATE  THERMAL  CONDUCTIVITY 

TCON=( (  0.002  66*(  CU)  **3.0  ) ) - ( 3 . 2* { C ( K ) **2 .0 ) )  +  ( 1073.33 
1*C ( K)  ) 

B+286000.0 ) /1000000.0 

C  CALCULATE  THE  PRANDTL  NO. 

PR= (VIS*3600.0) /TCON 

C  CALCULATE  THE  NUSSELT  NO. 

ANU=0.023*(RE**0 . 8 ) * ( PR**0 . 3 ) 

C  CALCULATE  THE  FILM  COEFFICIENT 

C( I )= (ANU*TCON*12.0)/(EQD  *3600.0) 

RETURN 

END 
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cm ue  3k i ruosaus 

(G#SI*SlV)\(A»£o* 
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).u-00Ct>  I\(O.0~  d8S  +  8 
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C/3 


APPENDIX  I-B 


Listing  of  CSMP  Model  Configuration 
and  Specification  Statements 


X 


A 


CUNT  I  NUOUS 
DIGITAL  ANALOG  S 


SYSTE  1  10  DELING 
MULATJR  PROGRAM 


PROGRAM 
FOR  THE 


IBM  1300 


REEL)  HELP  -  ENTER  1 


OTHERWISE  HIT  EOF 


IF  CONSOLE  DATA 
IF  POC  DATA  SWL 


OUTPUT  NAME 


SWITCH 

ES  TO  BE  US 

El) ,  ENTER 

0 

'CrlES  A 

RE  TO  BE  US 

ED,  ENTER 

1 

CON F IGOR A 

T 1  ON  SPECIFICATION 

.OCR 

TYPE 

INPUT  1 

INPUT  2 

1  NPUT 

1 

K 

0 

0 

0 

2 

1 

0 

10 

11 

3 

1 

0 

11 

12 

4 

1 

17 

12 

13 

5 

1 

0 

13 

0 

0 

+ 

-1 

2 

0 

7 

+ 

-2 

“T 

J 

J 

3 

+ 

-  j 

4 

0 

0 

+ 

4 

-5 

0 

10 

X 

68 

G 

0 

11 

»/ 

A 

14 

7 

0 

12 

V 

A 

15 

8 

0 

13 

>/ 

A 

1G 

0 

0 

14 

G  8 

2 

0 

13 

4 

GO 

4 

0 

1G 

5 

GO 

4 

0 

17 

W 

4 

20 

0 

13 

1 

U 

2  G 

27 

10 

1 

0 

2  7 

28 

20 

1 

•t  “T 

28 

29 

21 

1 

U 

20 

0 

22 

+ 

-2 

18 

0 

23 

+ 

-18 

10 

0 

24 

+ 

-10 

20 

0 

25 

+ 

20 

-21 

0 

2u 

y 

A 

63 

22 

0 

2  7 

V 

A 

3  0 

23 

0 

23 

X 

31 

2  4 

0 

20 

>/ 

A 

32 

25 

0 

30 

T 

D 

63 

18 

0 

31 

4 

G  9 

20 

0 

32 

5 

GO 

20 

0 

*»  •* 

D  J 

W 

20 

36 

0 

34 

1 

0 

42 

43 

3  5 

1 

0 

43 

44 

36 

1 

40 

44 

45 

37 

1 

0 

45 

0 

• 

IC/PAk  NAME 


56 

+ 

-13 

34 

0 

39 

+ 

-34 

35 

0 

40 

+ 

-35 

36 

0 

41 

+ 

36 

-37 

0 

4  2 

V 

A 

6  o 

33 

0 

43 

A 

46 

39 

0 

44 

X 

47 

40 

0 

45 

v 

A 

43 

41 

0 

46 

“* 

;> 

68 

34 

0 

4  7 

4 

GO 

30 

0 

43 

5 

69 

36 

0 

49 

U 

36 

53 

0 

50 

K 

0 

0 

0 

51 

1 

0 

59 

60 

52 

1 

0 

60 

61 

53 

1 

66 

61 

62 

54 

1 

0 

62 

0 

55 

+ 

-34 

51 

0 

56 

+ 

-51 

52 

0 

57 

+ 

-52 

53 

0 

53 

+ 

5;> 

-54 

0 

59 

/\ 

6  3 

55 

0 

60 

\/ 

/\ 

63 

56 

0 

61 

X 

6  4 

57 

0 

62 

X 

65 

58 

0 

63 

63 

51 

0 

6  4 

4 

69 

53 

0 

6  5 

5 

6  9 

53 

0 

6  6 

A 

53 

67 

0 

67 

K 

0 

0 

0 

6  3 

2 

50 

0 

0 

6  J 

K 

0 

0 

0 

INITIAL  CON 

:)  i  t  i  o:ms 

AMD  PARAMETERS 

BUCK 

1  C/PARI 

PAR  2 

PAR  3 

1 

57.3000 

0.0900 

0.0000 

2 

64.4000 

■0 . 400  0 

0.9435 

*T 

:> 

64. 4000 

■7.514 4 

8. 0673 

4 

147. 2000 

•0.4182 

-0.7363 

5 

147.2000 

5.4161 

0. 0000 

14 

0.3150 

0 . 0000 

0. 0000 

15 

1.5410 

0  .  8  7  5  0 

0.0000 

16 

1.5410 

0.8750 

0. 0000 

17 

-2. 0300 

2 . 0300 

0  .  0000 

13 

64.4000 

0.  40 0  0 

0.9435 

19 

64.4000 

- 

7.5144 

3  .  0  6  7  3 

20 

147. 2000 

- 

0.4182 

-0.7363 

21 

147. 2000 

5.4161 

0 . 0000 

30 

0. 3150 

0 . 0000 

0 . 0000 

31 

1.5410 

0. 8750 

0.0000 

32 

1.5410 

0.3750 

0 . 0000 

—9  ** 

J>  J 

-2. 0300 

2 . 0300 

0  ,  0  0  0  0 

3  4 

6  4 . 4  0  0  0 

- 

0 . 4  0  0  0 

0.9435 

5  4 .  40  0  J 

147.2000 
147. 2u0u 
0.8150 

1.5410 

1.5410 
-2.0500 
6.1500 
5  4 .  4  0  0  0 
64.4000 

147.2000 

147.2000 
0.8150 

1.5410 

1.5410 
-2.0500 
140. 8000 
15.0000 
5.0750 


-7.5144 
-0.4132 
5.4151 
0 . 0000 
0.8750 
0.8750 
2.0500 
0 . 0000 
-0. 4000 
-7.5144 
-0.4132 
5.4161 
0.0000 
0.8750 
0.8750 
2.0300 
0 . 0000 
20.0000 
0.0000 


3.0673 
-0.7353 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0 . 0000 
0.0435 
3.0673 
-0. 7363 
0  .  0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0 . 6500 
0.0000 
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APPENDIX  II-A 


Error  in  the  Trapezoidal  Rule  Approximation 

to  an  Integral 

Consider  integrating  the  function  f(t)  over  the  interval  [a,b] 
using  the  numberical  trapezoidal  approximation  method.  The  approxi¬ 
mation  error  may  be  shown  (10)  to  be 

E  =  (b-a)  At2  f"  (n) ,  a  <  n  <  b  ...  II-A-1 

12 

Assume  a  time  record  length  of  25  time  units,  500  data  points,  and 
f"  (n)  equals  unity,  then 

E  =  0 (At2)  =  25_  (0.05)2 

12 

which  is  in  the  order  of  a  1  per  cent  error.  This,  of  course,  would 
be  quite  acceptable  for  most  engineering  applications.  The  method 
breaks  down,  however,  because  the  quadrature  products  to  be  evaluated 
are  composed  of  sine  or  cosine  terms,  which  oscillate  faster  as 
frequency  increases.  This  results  in  the  rate  of  change  of  the  slope 

evaluated  at  n,  i.e.,  f"  (ri)  becoming  much  larger  than  unity  as  the 
oscillation  rate  increases.  This  gives  rise  to  a  corresponding 
increase  in  the  error  of  the  approximation. 


' 
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APPENDIX  II-B 


Filon's  Method 


To  preserve  the  accuracy  of  the  Filon  formulae  when  0  is  less 
than  0.35  radians,  the  following  equations  should  be  used  for 
evaluating  a,  3,  and  y: 

a  =  26^  -  20i  -  l£_  -  .  ...  u-b-i 

45  315  4725 

3  =  _2  +  262  -  40^  +  20®  -  40®  + .  ...  H-B-2 

3  15  105  567  22275 

( 

y  =  4  -  2§£  +  0^_  -  06  +  0 8  -  .  ...  ii-b-3 

3  15  270  11340  997920 


These  formulae  prevent  loss  of  significant  digits  and  resulting 
roundoff  errors. 


■ 
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APPENDIX  II-C 


FFT  Algorithm 

The  power  of  the  FFT  algorithm  is  due  to  the  possibility  of 
employing  a  formula  which  computes  a  transform  of  N  data  points, 
given  the  transform  of  N/2  data  points,  without  quadrupling  the 
effort;  one  would  expect  a  quadrupling  effect  from  the  considera¬ 
tion  that  N2  operations  are  required  by  an  N  point  transform  (21). 
The  formula  employed  is 


Z 


k 


\  +  Bk 


. . .  II-C-1 


Zk+N/2  Ak  Bk  W 


...  II-C-2 


where  k  =  0,  1,  2, 


,  N  -  1 
2 


and  W  =  exp  -j  1_ 

N 


...  II-C-3 


N/2  -1 


Z  f2i  (W  } 
i  =  0 


2v  ik 


. . .  II-C-4 


B. 


N/2  -1 

Z  f2i+l^W<^ 
i  =  0 


2x  ik 


. . .  II-C-5 


where  k  =  0,  1,  2, 


,  N  -  1 
2 


As  the  formulae  indicate,  an  N  point  time  history  is  considered 
in  two  parts,  the  even  indexed  values  (f_.)  and  the  odd  indexed 

2i 

values  (^21  +  1^*  individual  transforms  of  these  N/2  point 

sequences  are  computed  and  then  combined  by  using  Equations  II-C-1 
and  II-C-2.  This  basically  requires  an  additional  N  multiply-add 
operations.  Suppose  one  now  divides  N/2  sequences  to  obtain  four 
N/4  point  sequences;  this  would  require  an  extra  N/2  multiply-add 
operations.  Carrying  this  procedure  to  the  limit  until  arriving  at 
a  2-point  transform  beginning  with  2n  data  points,  one  will  cascade 
through  n  stages,  each  requiring  N  operations,  giving  a  total  of 
nN  operations  in  all.  This  compares  with  N2  operations  by  a 
straightforward  procedure. 
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APPENDIX  III-A 


Listing  of  Process  Programs 


3. A 


1 


C 

C 

c 

c 

c 

c 

c 

c 

c 


********  ***********##*#.*.#.*  **#**#***#*##************.£#.£ 


PROCESS  PROGRAM  HTBAL 


* 


* 


EXTERNAL  GENER 

DIMENSION  TFLO( 100) ♦ SFLO ( 100)  *  T I T  <  100)  *TOT(100) 
DIMENSION  SIT(IOO) *SOT( 100) 

COMMON  1111(700) 

COMMON  LUR*LUW*LUP*NDATA*NSS*NPUL  * IPULW*DT 
COMMON  Y ( 400 ) 

C  DEFINE  LOGICAL  UNIT  NUMBERS 


LUR=  1 
LUW  =  2 
LUP  =  6 


C 


DEFINE  LOOP  I • D* •  S 


I DTF=  38  4 
I DSF=  38  5 
IDT  I T  =  386 
I DTOT=387 
I DS I T  =  3  88 
I DSOT  =  3  89 

C  ***  KEYBOARD  REQUEST  *** 

WRITE  ( LUW  *  5 ) 

5  FORMAT  ( 1  OX  * ' TYPE  I  SEC  *  ) 

CALL  FFINP (LUR* 1 *0# ISEC* IEROR ) 

C  CHECK  SIZE  OF  ISEC 

IF  ( I  SEC- 100 )  6*6*35 

6  CONTINUE 
TF=0 • 0 
SF=0 • 0 

T I =0  •  0 
T0  =  0  •  0 
SI=0.0 
S0=0.0 


.. 


L)  -IT  .  1  ..  I 

, 

j,  W  J» 


Ji . — ^  1  J  •  •  .  1  ' 


♦  .a.  i  iooj  :  -  :  bb a 


4i8£=*3T0I 

V8£*TOTOI 
'  *  t  =  T  I  • '  1 
^8£=TOaOI 


(£«WUJ)  3  T  1 5R  a 

(  »  3  ;t:  i  qy  *  *  ■  OX  )  t  t 

(  ,-.3 :  *  o  u*  •  U  -i )  ^  J A  5 


oiei  bo  3iie  >o3ho 
♦  a.  .  .  0  >1-0321  } 


. 

0.0=  n 

. 

. 


3. A 


2 


C  INITIALIZE  TIME 

CALL  T IME ( KHOUR *KMIN»KSEC ) 

KT  I ME-KM I N*60+KSEC 
DO  20  1  =  1  *  I  SEC 

10  CALL  T I  ME ( JHOUR * jM I N * JSEC ) 

JT I ME= ( ( JHOUR -KHOUR )*360Q)+(JMIN*60)+JSEC 
NDIFF=KT IME+I-jT IME 
IF  (NDIFF)  15*15*10 

C  TIME  TO  READ  THE  VARIABLES  AGAIN 

15  CALL  GTVLU ( IDTF* 1 *TFLO( I )  *  IERR # IF ) 

CALL  GTVLU ( IDSF*1.SFLQ( I ) * IERR* IF) 

CALL  GTVLU ( I  DTI T  *1 ,TIT ( I )  *  IERR* IT ) 

CALL  GT VLU ( IDT0T*1*T0T(I ) *IERR*IT) 

CALL  GTVLU( IDSIT  *1*SIT (  I  )  .  IERR, I T ) 

CALL  GTVLU ( IDSOT  *  1 ,SOT ( I  )  , I  ERR , I T ) 

20  CONTINUE 

C  CALCULATE  THE  AVERAGES 

DO  22  1  =  1*1  SEC 
TF=TF+TFLO ( I ) 

SF=SF+SFLO( I ) 

T I =  T I +T I T ( I ) 

TO=TO+TOT ( I ) 

S I =S I +S I T ( I ) 

SO=SO+SQT ( I ) 

22  CONTINUE 
TF  =  TF / I  SEC 
SF  =  SF / 1  SEC 
T I =T I / I  SEC 
TO  =  TO/ I  SEC 
S I =S I / I  SEC 
SO  =  SO/ I  SEC 

C  CALCULATE  THE  HEAT  BALANCE 

RAT  10= (SF* (SI -SO )  ) / ( TF*( TO-TI  )  ) 

C  WRITE  OUT  THE  RESULTS 

WRITE! LUW  *  25 ) 

25  FORMAT (////// *20X» 'HEAT  BALANCE  DATA') 

WRITE  (LUW*30)  TF*TI*TO 
30  FORMAT  (///*'TUBE  FLOW '  * F5 . 2  ♦  '  U . S . G • P . M . '  » 3X  , 
l'INLET  TEMP'  *F5. 1  * 'DEG*F 1  *3X» 'OUTLET  TEMP'* 
2F5.1*  ' DEG • F '  ) 


•>  • 


* 


<  >>•  !I  vr  » S  J  i  H'i  J  1/ 

..  ‘  -  "  I  v  : 

; « I  I  -  o< 

(  D  U*  I  u*  ^’JCK, )  j  I  T  JJ/O 
L )  ’  •  •••  -  1,1  C>  =  "  M  !  L 

I  tl  - i  •  ‘IT  ill  :/■ 

i .  ej  *  i  mu  ) 


;v  I A  .  J>  1 9  W  '•  T  OA  :  (•  1  r  ’  I  T 

:  . ,  :  ,{j.  r .  . .  •  )'u  m 

;  . ,  •  (  i  •  .  : 

•  ,  ...  ITOI 1UJVT3  J  IA 

:  I  *  !  1  )  T  :T  »  ;  ♦  T  T  a  I  )  J  J  V T £)  J  -  AD 

(  I  I  .  'ii*  m  T  !  .*  X  *  T  MOD  .OVT  JJ/O 

(  ;  .  t  (  i  )  .  .  *  T(  i  JVTc  J  J  ■  D 

^  I  V-  D 


£3cDAfl  3  V  A  NT  3TAJJDJAD 

:  N 

:  ;  •  i+-T=  -: 

(  1  )0J  *•"  • 

(  I  )T1T+'1T*IT 
t I 1 TOT+OT-OT 

1 1  )Tu+ia*u 

D3eiM8*12 
.  :  i  .  : :  r 
D32 I \0T  *0T 


i 


i*  ’  j  ‘  ♦  \\  \  )  T  A 


i  »  .  •  '  .  '  «  I  M 


,  •  r  J  *  : 


n  n 


3. A  3 


WRITE  ( LUW  *  3  5 )  SF»SI*SO 

33  FORMAT  ( / / *  *  SHELL  FLOW »  » F5 . 2 » ' U . S . G . P . M . •  t 3X » 
1 ' INLET  TEMP*  *F5. 1  * *DEG.F '  *3X* 'OUTLET  TEMP'» 
2F5 • 1 » 1  DEG • F 1  ) 

WRITE  ( LUW  *  40 )  RATIO 
40  FORMAT  ( / / / *  1  OX  *  ' QS /QT = '  * F6 • 3 ) 

CHECK  DATA  SWITCHES  -  ABORT  WITH  5*CHAIN  WITH 
SIX  TO  PROGRAM  GENER 

CALL  DATSW(5* IPOS) 

IF  (IPOS-1)  5  5  >  5  5  »  45 
45  CALL  DATSW(6 *  IPOS) 

IF  (IPOS-1)  50  *  50  »  6 
50  CALL  CHAIN (GENER ) 

55  CALL  VIAQ 
END 


•  ' 


•  » i  (t  *  J)  ; r 

•  ’  * •  i 

.  •  .  O'*  .  ■ ; 

^  •  !  )  1 

•  1  1  !  •  .  *  ■ 


t;  ;  AT  AC  *3:’  O 

{  ;  I •  )  ' /  J -  A3 

.  »  C  »  .  *  *  i  >  3  1 

. 

■  .  c  •  c  i  -  I  )  I 


10 

15 

20 

25 

30 

35 

100 


3  •  A 


4 


#  *  *  ***********  *  **********  *  *  *  **********  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 


* 

* 

*  PROCESS  PROGRAM  GENER 

*  - - - - 

* 

* 


*  *  *  *********************  *  ***********  *  **********  * 


* 

* 

* 

* 

* 

* 

****** 


EXTERNAL  TPLOT 
COMMON  1111(700) 

COMMON  LUR»LUW»LUP.NDATA»NSS*NPUL, IPULWfDT 
COMMON  Y ( 400 ) 

***  KEYBOARD  REQUEST  *** 

WRITE  ( LUW  » 1 0 ) 

FORMAT (//// *  *  TYPE- I  SUB  * IPULW* I  RATE •  ) 

CALL  FFINP  ( LUR  *  3  *  0  *  I  SUB  *  0  *  I PULW  »  0 ♦ I RA TE  *  I EROR ) 

CALC  FORCING  FUNCTION  PARAMETERS 

NDATA=40* IRATE 
NSS=5* IRATE 
NPUL= IPULW*IRATE 
DT=1. 0/IRATE 

CALL  APPROPRIATE  GENERATING  SUBROUTINE 

GO  TO  ( 15 >20 *25 >30>35 ) * ISUB 

CALL  PULI 

GO  TO  100 

CALL  PUL 2 

GO  TO  100 

CALL  PUL 3 

GO  TO  100 

CALL  PUL4 

GO  TO  100 

CALL  PUL 5 

CALL  CHAIN ( TPLOT  ) 

END 


_ _ _ _ _ 
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SUBROUTINE  PULI 


COMMON  1111(700) 

COMMON  LUR»LUW»LUP*MDATA»NSS*NPUL* IPULW*DT 
COMMON  Y ( 400 ) 

C  ***  KEYBOARD-PULSE  IDENTIFICATION  *** 

WRITE  ( LUW  *10) 

10  FORMAT  ( / / // * 10X  * •  RECTANGULAR  PULSE*) 

C  GENERATE  THE  FORCING  FUNCTION 

DO  20  I  =  1 » NDATA 

Y  {  I  )  =0.0 
20  CONTINUE 

IBGN=NSS+1 
I  END= I BGN+NPUL 
DO  25  I  =  I  BGN  » I  END 

Y  (  I  )  =  1  •  0 
25  CONTINUE 

RETURN 

END 
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SUBROUTINE  PUL2 


COMMON  1111(700) 

COMMON  LUR»LUW*LUP*NDATA>NSS*NPUL* IPULW*DT 
COMMON  Y ( 400 ) 

C  ***  KEYBOARD-PULSE  IDENTIFICATION  *** 

WRITE  ( LUW  *10) 

10  FORMAT  (////* 10X *' RAMP  PULSE') 

C  GENERATE  THE  FORCING  FUNCTION 

DO  20  I  =  1  • NDATA 
Y  (  I  )  =0.0 
20  CONTINUE 

I  BGN  =  NSS+ 1 
I  END= I BGN+NPUL 
DELY= 1 • 0/NPUL 
SUM=0 . 0 

DO  25  I  =  I BGN  *  I  END 
Y ( I ) = Y ( I ) +SUM 
SUM=SUM+DELY 
25  CONTINUE 
RETURN 
END 
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SUBROUTINE  PUL3 


COMMON  1111(700) 

COMMON  LUR*LUW*LUP»NDATA*NSS*NPUL  * IPULW*DT 
COMMON  Y ( 400 ) 

N.B.  AN  EVEN  NO.  OF  SECONDS  SHOULD  BE  SPECIFIED 
FOR  THIS  PULSE 

***  KEYBOARD-PULSE  I  DENT  I F I  CAT  I OC  *** 

WRITE  ( L  U  W  » 1 0 ) 

10  FORMAT (////*10X* ' TRIANGULAR  PULSE') 

C  GENERATE  THE  FORCING  FUNCTION 

DO  20  1=1 *N DATA 
Y  (  I  )  =0.0 
20  CONTINUE 

I  8GN=NSS+ 1 
I  END= I BGN  +  NPUL 
DELY=2 • 0/NPUL 
SUM=  0.0 

I M I D  = I 8GN  + (  ( I  END- 1 BGN ) / 2 ) 

DO  35  I  =  I  BGN  *  I  END 
Y ( I ) = Y ( I ) +SUM 
W  Y  =  Y  (  I  ) 

IF  (I-IMID)  25*30*30 
25  SUM=SUM+DELY 
GO  TO  35 

3  C  SUM=  SUM~D  ELY 
35  CONTINUE 
RETURN 
END 
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SUBROUTINE  PULA 


COMMON  1111(700) 

COMMON  LUR*LIJW*LUP>NDATA»NSS»NPUL*IPULW*DT 
COMMON  Y ( 400 ) 

C  ***  KEYBOARD-PULSE  IDENTIFICATION*** 

WRITE  (LUW*10) 

1C  FORMAT  (////* 10X »' DISPLACED  COSINE  PULSE’) 

C  GENERATE  THE  FORCING  FUNCTION 

DO  20  I  =  1 »  NDATA 
Y (  I  )  =0.0 
20  CONTINUE 

I BGN=NSS+ 1 
I  END= I BGN  +  NPUL 
L  =  0 

DO  25  I  =  I BGN  *  I  END 

Y  (  I )= (1.0-COS(6.2  8*L*DT/IPULW)  )/2.0 
L  =  L+ 1 

25  CONTINUE 
RETURN 
END 
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SUBROUTINE  PUL5 


COMMON  1111(700) 

COMMON  LUR»LUW*LUP*NDATA*NSS»NPUL » IPULW*DT 
COMMON  Y ( 400 ) 

C  ***  KEYBOARD-PULSE  IDENTIFICATION  *** 

WRITE  ( LUW  *10) 

10  FORMAT  (////♦ 10X »' RECTANGULAR  DOUBLET  PULSE*) 

C  GENERATE  THE  FORCING  FUNCTION 

DO  20  I = 1 *NDATA 
Y (  I  )=Q.Q 
20  CONTINUE 

I BGN=NSS+ 1 
I END= I BGN+NPUL/2 
DO  25  I  =  I BGN  *  I  END 
Y  (  I  )  =  1  •  0 
25  CONTINUE 

I BGN= I END+1 
IEND=NSS+1+NPUL 
DO  30  I  =  I BGN  *  I  END 
Y ( I ) =-1.0 
30  CONTINUE 

RETURN 
END 
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C  ****************************************************** 


c  *  * 

c  *  * 

C  *  MAINLINE  PROGRAM  TPLOT  * 

C  *  * 

C  *  * 

C  *  * 


C  *************  ******  **************  ***********  ********** 

c 

EXTERNAL  PULSE 
DIMENSION  X ( 400 ) 

COMMON  1111(700) 

COMMON  LUR  * LUW  *  LUP ♦ NDAT A  *  NSS  » NPUL  » IPULW»DT 
COMMON  Y ( 400 ) 

CALL  PLOTDd) 

5  I DEV=7 

C  GENERATE  THE  X  VECTOR 

SUM= 1 • 0 

DO  10  I  =  1  *  NDATA 
X { I ) =SUM 
SUM=SUM+1 .0 
10  CONTINUE 

C  DEFINE  SCALE  FACTOR  &  STARTING  PEN  LOCATION 

XSF=8.Q/NDATA 
YSF=5. 0/4.0 
XPPL=0 • 0 
YPPL=-4 .0 

CALL  SCALF(XSF*YSF*XPPL»YPPL) 

DEFINE  THE  PLOTTING  GRID 
***  X  AXIS  *** 

RL  =  NDAT  A/  16 • 0 

CALL  FGRID(0*0.0»-2.0*RL*16) 

***  Y  AXIS  *** 

CALL  FGRID(1*0.0*“2.0*0.5*8) 

DEFINE  THE  DESIRED  ANNOTATION 
***  X  AXIS  *** 


. 
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XN=NDATA*0. 2/8*0 

CALL  FCHAR  ( -XN * -2 . 2 5  *  0 . 2  *0.2 *0.0  } 

WRITE  (  I  DEV  *  2  0 ) 

20  FORMAT  (•  0 '  * 8X *  1 1 0 '  . 8 X  * ' 20 '  *  8 X  *  *  30 1  * 8 X  * ' 40 '  ) 

XN=NDATA*2.35/8. 0 

CALL  FCHAR (XN *-2 .6 >0.3 *0.2 *0.0 ) 

WRITE  ( I  DEV  *  30 ) 

30  FORMAT (  • TIME  (SECS )  '  ) 

C  ***  y  AXIS  *** 

XN=NDATA*0. 20/8.0 
YN  =  2 • 2 

CALL  FCHAR ( -XN*-YN*0.2 *0.2  ♦  1 . 57 ) 

WRITE  ( I  DEV  *  40 ) 

40  FORMAT  ( ■ -2 ' *4X * ' -1 • *4X * •  0 ' * 4X » • +1 ' *4X » ' +2 1 ) 

C  NOW  PLOT  THE  DATA 

CALL  FPLOT ( +1 »0. 0*0.0 ) 

CALL  FPLOT (-2*0. 0*0.0) 

DO  50  I = 1 *NDATA 
X  N  =  X  (  I  ) 

YN  =  Y (  I  ) 

CALL  FPLOT (0 *XN*YN) 

50  CONTINUE 

CALL  FPLOT ( +1 *XN *-4.0 ) 

C  ***KE YBOARD  REQUEST*** 

WRITE  ( LUW  *  60 ) 

60  FORMAT  (10X*»IF  PLOT  DESIRED  TYPE  A  1»» 

1  •  -  IF  NOT  TYPE  A  2 '  ) 

CALL  FFINP (LUR* 1 *0* IFLAG* IEROR) 

IF  (  I  FLAG- 1 )  70  *70*80 
70  CALL  PLOTD(O) 

GO  TO  5 

80  CALL  CHAIN (PULSE  ) 

END 
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C  ******#****##*#**#*##*#****#•■********•*•«■**#■■*#■**#■***##*** 


C  *  * 

c  *  * 

C  *  PROCESS  PROGRAM  PULSE  * 

C  *  * 

c  *  * 

c  *  * 


C  *********##*#*****#***•***#•******#*  #  -;*•  *  -x-  **************  *  * 

C 

EXTERNAL  ASORT 
COMMON  1111(700) 

COMMON  LUR  > LUW * LUP > NDATA > NSS » NPUL f IPULW»DT 
COMMON  Y ( 400 ) 

C  ***  KEYBOARD  REQUEST  *** 

WRITE  ( LUW  *10) 

10  FORMAT  (////*' TYPE-SPAN  (OF  PULSE  IN  PERCENT)’) 

CALL  FFINP  (LUR*1>1*SPAN» IEROR) 

C  DETERMINE  PRESENT  VALUE  OF  OUTPUT  TO  THE  VALVE 

IDTF=384 

CALL  GTVLU  (  I DTF *  3 » YBASE » I  ERR > JDESC ) 

C  WRITE  OUT  PRESENT  VALUE 

WRITE  ( LUW  >  20 )  YBASE 

20  FORMAT  (’PRESENT  OUTPUT  TO  THE  VALVE  ='.F6.2*'(') 

C  CALCULATE  THE  Y  VALUES  IN  PERCENT 

DO  30  I  =  1 »  NDATA 
Y (  I  )=YBASE+(SPAN*Y( I )  ) 

30  CONTINUE 

NBGN=NDATA+1 
NF  I  N=NDATA+50 
DO  35  I =NBGN  »  NF I N 
Y  (  I  )  = YBASE 

35  CONTINUE 

WRITE  ( LUW  >  36 )  NDATA 

36  FORMAT  (’NO.  OF  PULSING  POINTS  = 1 *16) 

C  PLACE  TUBE  FLOW  CONTROL  LOOP  ON  MANUAL 

CALL  NONOP (0180) 

CALL  KILL 
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C  MAKE  VCDAQ  SUBROUTINE  IDLE 

I  FLAG  =  0 

CALL  VCDAQ ( I  FLAG ) 

C  START  IMPLEMENTATION  OF  THE  PULSE 

CALL  CYCLE ( 9  >  7  *  2 ) 

C  START  HIGH  SPEED  DATA  ACQUISITION 

I  FLAG  =  2 

KEEP  CHECKING  TO  SEE  WHETHER  THE  40  SECOND 
DATA  ACQUISITION  PERIOD  IS  FINISHED 

40  CALL  TVCDQ (  ICHEK) 

GO  TO  (40>50)  > ICHEK 
50  WRITE  ( LUW  *  60 ) 

60  FORMAT  ('HIGH  SPEED  DATA  ACQUISITION  COMPLETED') 
CALL  CANCL ( 7 ) 

C  PLACE  TUBE  FLOW  CONTROL  LOOP  BACK  ON  AUTOMATIC 

CALL  OPER ( 0 1 80 ) 

CALL  CHAIN(ASORT) 

END 
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SUBROUTINE  VALVE 


THIS  SUBROUTINE  ACCESSES  THE  DAC  ,  PROVIDES  IT  A 
Y  VALUE, THEN  THE  TUBE  FLOW  CURRENT  OUTPUT 
STATION  IS  PULSED. THIS  SUBPROGRAM  IS  CALLED  5 
TIMES  PER  SECOND 

DIMENSION  I  DATA ( 3  )  » I PO ( 3 ) 

COMMON  1111(700) 

COMMON  LUR > LUW * LUP , NDAT A »NSS ,NPUL » IPULW*DT 
COMMON  Y ( 400 ) 

DATA  11/0/ 

11=11+1 

I0UT=32767.0*Y( I  I  )  / 100.0 

C  ACCESS  THE  DAC 

IDATAtl )  =  I  OUT 
I  DATA ( 2  )  =  1 
IDATA ( 3  )  =  2 

5  CALL  DAC ( 10001  , I  TEST ) 

GO  TO  (5,10), ITEST 

10  CALL  DAC ( 11001, I  DATA (1)  ,IDATA(3)  ) 

I PO (  1  ) =512 
I  PO ( 2 )  =  120 
I PO ( 3 ) =2 

CALL  PO (  12001  , IPO ( 1 )  ,IPO(3)  ) 

CALL  INTEX 

END 
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C 

C 

C 

C 

C 

c 

c 

c 

c 


*  *  *  *  *  *  -x-  *  -X-  *  *  -X-  *  *  *  *  *  *  -X  -X-  *  -X-  -if  -x-  *  -x-  *  *  -x-  ******  *  *  *  *  *  *  *  *  *  *  *  x-  -x-  *  *  *  *  *  x 


* 

* 

* 

* 

* 


PROCESS  PROGRAM  ASORT 


* 

* 

* 

* 

* 

* 


*  -x-  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *****  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 


INTEGER  OUT FL (83 ) 

EXTERNAL  VTFLO 

DIMENSION  IFUNC ( 8 )* INFLE ( 323 ) 

COMMON  I  I  I  I  ( 700 ) 

COMMON  LUR  *  LUW  *  LUP  *  NDAT  A  *  NSS  *  NPUL  * IPULW»DT 
COMMON  Y ( 400 ) 

DEFINE  FILE  5  1  ( 40 » 80 > U » JO )  » 52  (  40  > 80 » U » JO ) 
DEFINE  FILE  65 ( 4 1 » 320  ♦ U *  JO ) 


C  DEFINE  NECESSARY  PARAMETERS  FOR  SUBROUTINE  FLSRT 


C  SORT  THE  FLOW  DATA 

IFUNC ( 1 ) =  040  9  6 
IFUNC (2 ) =40 
INFLE(1 ) =  6  5 
INFLE (2 ) =320 
INFLE ( 3 ) =41 
OUTFL ( 1 ) =  5 1 
OUTFL ( 2 ) =80 
OUTFL ( 3 ) =40 

WRITE  (LUW >10) 

10  FORMAT  ( / / / ♦ 5  X  » 'NOW  SORTING  FLOW  DATA') 
CALL  FLSRT ( IFUNC » INFLE »QUTFL ) 

C  WRITE  OUT  ERROR  PARAMETERS 

WRITE  ( LUW  >  2 0 )  (  I FUNC ( J )  ♦ J  =  1 » 8 ) 

20  FORMAT  ('ERROR  PARAMETERS= '  * 8  I  5 ) 

C  SORT  THE  TEMPERATURE  DATA 

IFUNC ( 1 ) =04098 
OUTFL ( 1 ) =52 
WRITE  (LUW *30) 


—  —  - - — - - 
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30  FORMAT  (///»5X»»NOW  SORTING  TEMP 
CALL  FLSRT ( IFUNC ♦ INFLE  »OUTFL ) 

C  WRITE  OUT  THE  ERROR  PARAMETERS 

WRITE  ( LUW  *  20 )  I FUNC ( 3 )  ♦  I  FUNC  l  4 ) 
CALL  CHAIN(VTFLO) 

END 


DATA '  ) 
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*******************  **-B-***tf*##*tf-8-#*#*******##****  ****** 


*  * 

*  * 

*  PROCESS  PROGRAM  VTFLO  * 

*  - - - — - —  * 

*  * 

*  * 


*  -a-#-****-*-*  **********  *******************  **************** 


EXTERNAL  MVTMP 

DEFINE  FILE  5 1 ( 40  * 80 »U * JO )  » 52 ( 40  * 80  * U # JO ) 

DEFINE  FILE  53 ( 40  *  80 #U * JO ) 

DIMENSION  IVOLT ( 80) #FLOW ( 80) 

COMMON  II  I  I  (700) 

COMMON  LUR*LUW*LUPtNDATA*NSS*NPUL*IPULW*DT 
COMMON  V ( 400 ) 

JREC- 1 

DO  30  L= 1  * 20 

READ  ONE  RECORD  OF  FLOW  DATA  (I.E.  80  INTEGER  VALUES) 

READ  ( 5 1  *  L )  IVOLT 
DO  10  I  *  1 »80 

CONVERT  INTEGER  VALUES  TO  REAL  VALUES 
V=FLOAT ( IVOLT ( I ) ) 

CONVERT  TO  ACTUAL  VOLTAGE 
V=5.0/32767«0*V 

CONVERT  TO  ACTUAL  FLOW  (3/4  INCH  TURBINE  ASSUMED) 

FLOW ( I )=( ( (0«916809*V~ 13*2838 )*V+287*Q47 )*V+11»0247 ) 
1*0*0195 
10  CONTINUE 

WRITE  THE  TWO  RECORDS  OF  REAL  VALUES  TEMPORARILY 
IN  FILE  53 

KBGN=1 
KEND=40 
DO  20  J* 1  * 2 

WRITE  (53* JREC)  (FLOW(K) » K=KBGN ♦ KEND ) 

KBGN=KBGN+40 
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KEND=KEND+40 
JREC= JREC+1 
20  CONTINUE 
30  CONTINUE 

C  NOW  TRANSFER  VALUES  FROM  FILE  53  TO  FILE  51 

DO  40  JRECS 1 *40 

READ  (53‘JREC)  ( FLOW < J ) * J= 1 »40 ) 

WRITE  ( 5 1 ' JREC )  ( FLOW ( J ) * J* 1 *40 ) 

40  CONTINUE 

WRITE  ( LUW  *  50 ) 

50  FORMAT  ( / // *  5X  *  * FLOW  VALUES  CONVERTED  TO  ENG.  UNITS1 ) 

C  CHECK  FIRST  RECORD  OF  DATA 

READ  <  5 1  '  1  )  ( FLOW ( J ) * J  =  1  *  40 ) 

SUM=0 . 0 
DO  60  J* 1 *40 
SUM=SUM+FLQW ( J ) 

60  CONTINUE 

SUM=SUM/40 .0 
WRITE  ( LUW  *  70 )  SUM 

70  FORMAT  (*AVG.  FLOW  FOR  1  ST  SEC  OF  RUN=»*F5.2) 

WRITE  (LUP*65) (FLOW(J) »J=1*40) 

65  FORMAT ( 10F7 . 3 ) 

CALL  CHAIN(MVTMP) 

END 
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**#****#*#**#*#***#-******##*###***■#****■»•■**#**■***#**#** 


*  * 

#  * 

*  PROCESS  PROGRAM  MVTMP  * 

*  - - — * 

*  * 

*  * 


************  *******  *•**#*******■«■  *  *********  **  ****  *  *****  * 


EXTERNAL  BPLOT 

DEFINE  FILE  5 1 ( 40 * 80 *U . JO ) *  52 ( 40 * 80 *U * JO ) 

DEFINE  FILE  53 ( 40 *80 *U * JO ) 

DIMENSION  RBT ( 5 ) *REF<5) *TC(5) *MVOLT(80) * TEMP (80) 
COMMON  I II  1(700) 

COMMON  LUR*LUW*LUP*NDATA*NSS*NPUL*IPULW*DT 
COMMON  Y ( 400 ) 


C  DEFINE  PERTINENT  LOOPS 


I  DRBT  =  392 
I  DREF  =  393 
I DTC=400 


C  GET  AVG.  VALUES  OF  ABOVE  LOOPS  OVER  A  5  SEC.  PERIOD 

DO  30  1=1*5 

CALL  TIME ( JHOUR* JMIN* JSEC ) 

JSEC= JMI N*60+JSEC 
10  CALL  T IME ( I  HOUR  * IMIN*ISEC) 

ISEC=IMIN*60+ISEC 
IDIFF=ISEC-JSEC 
IF  (IDIFF-1)  10.20*20 
20  JSEC= JSEC+1 

CALL  GTVLU ( I DRBT *  1  * RBT ( I ) *IERR.IU) 

CALL  GTVLU( ID REF  *1*REF( I ) .IERR.IU) 

CALL  GTVLU ( IDTC*1*TC( I ) *IERR*IU) 

30  CONTINUE 

C  AVERAGE  THE  READINGS 

VRBT=0 • 0 
VREF=0 • 0 
VTC=0 • 0 
DO  40  1=1*5 
VRBT*VRBT+RBT( I ) 

VREF=VREF+REF( I ) 

VTC=VTC+TC( I ) 
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40  CONTINUE 

VRBT=VRBT/5.0 

VREF=VREF/5.0 

VTC=VTC/5.0 

DIVIDE  AVG  REF  VOLTAGE  BY  10  TO  ACCOUNT  FOR  THE 
DIFFERENCE  IN  EDIT  CHARACTERS  IN  THE  LOOPS 

VREF=VREF/10.0 

CALCULATE  THE  TEMP  IN  THE  COLD  JUNCTION  REFERENCE 
BOX  USING  A  SYSTEM  SUPPLIED  EQUATION 

TR=51.876*VRBT/VREF+41.Q 

CALCULATE  THE  COLD  JUNCTION  COMPENSATED  MV  FOR 
THE  CH-CO  T/C  IN  THE  REFERENCE  BATH 

VCOR= ( <0. 150726E-04*TR ) +0 • 3 1467 16E-0 1 ) *TR-1 *022466 

CORRECTED  T/C  MV  IS 

RFTMV-VCOR+VTC 

CONVERT  THE  T/C  MV  TO  AN  EQUIVALENT  TEMP  VALUE 
2=RFTMV 

IF  (2- 4.26298)  50*50*60 

50  T* ( t ( 0.83 5171 8E-02*2> -0.4 3445 16) *2+30. 8 3 3308) *2+ 
132*001089 
GO  TO  70 

60  T- ( ( ( 0.2 15048 lE-02*2) -0*2 575 562 1*2+29.656109 ) *2+ 
134.299718 

WRITE  OUT  THE  REFERENCE  BATH  TEMPERATURE 
70  WRITE  ( LUW  *  80  )  T 

80  FORMAT  (5X*‘TEMP  OF  REFERENCE  BATH  *»*F7.3) 

CONVERT  TEMP  OF  REF  BATH  TO  A  MV  EQUIVALENT 
FOR  A  CU-CO  T/C 

IF  (T-91.0)  90*90*100 

90  B= ( ( < 0.148 29 09E-07*T ) +0.108 1611 E-04 ) *T+0 . 2067943E-01 ) 
1*T-Q . 6737249 
GO  TO  130 

100  IF  (T-151.0)  110*110*120 

110  B=( ( ( 0.7703282E-08*T ) +0 . 9890805E-05 ) *T+0 • 2 1Q4733E-01 ) 
l*T-0. 6942869 
GO  TO  130 
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120  B  = ( ( (0.1301623E-07*T)+G.4450828E-05)*T+0.2235491E-01 ) 
l*T-0. 78613748 
130  CONTINUE 

C  CONVERT  THE  TEMP  DATA 

JREC= 1 

DO  200  L=  1  *  20 

READ  ONE  RECORD  OF  TEMP  DATA 
-I.E.  80  INTEGER  VALUES 

READ  ( 52 1 L )  MVOLT 
DO  180  1=1*80 

CONVERT  INTEGER  VALUES  TO  REAL  VALUES 
V  =  FLOAT ( MVOLT ( I )  ) 

CONVERT  TO  MV  (NOTE  5  PROBES  IN  PILE) 

X=V/32767. 0*20. 0/5.0 
ADD  REF  BATH  COMPENSATION 


X  =  X+B 


CONVERT  TO  TEMP  IDEG.F)  USING  CU-CO  T/C  TABLES 
IF  (X-1.309)  140*140*150 

140  TEMP (  I  )  =  (  (  <0.6025596E-02*X)-1.25243)*X+46.69154)*X 
1+32.01969 
GO  TO  180 

150  I F ( X-2 • 736  )  160*160*170 

160  TEMP  (  I  )  *=  (  (  (  0. 1 70  2  39  7E-0 1  *X  ) -1.04606  )*X+46. 05431  )  *X  + 
132.48098 
GO  TO  180 

170  TEMP ( I)«< ( ( -0. 65 9 3423 E-02*X) -0.6461 39 )*X+44. 35243 )*X 
1+34.63365 
180  CONTINUE 

WRITE  THE  TWO  RECORDS  OF  REAL  VALUES  TEMPORARILY 
IN  FILE  53 

KBGN= 1 

KEND=40 

DO  190  J*  1  ♦  2 

WRITE  (53' JREC)  ( TEMP ( K ) *  K  =  KBGN  *KEND ) 

KBGN=KBGN+40 

KEND=KEND+40 
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JREC= JREC+1 
190  CONTINUE 
200  CONTINUE 

C  NOW  TRANSFER  DATA  FROM  FILE53  INTO  FILE  52 

DO  210  JREC=1*40 

READ  ( 5 3  *  JREC )  ( TEMP ( J ) * J= 1 *40 ) 

WRITE  (52'JREC)  ( TEMP ( J ) * J= 1 ♦ 40 ) 

210  CONTINUE 

WRITE  ( LUW  *  220 ) 

220  FORMAT  ( /// *  5X  *  1  TEMP  DATA  CONVERTED  TO  ENG.  UNITS1) 

C  CHECK  FIRST  RECORD  OF  DATA 

READ  (52*1)  ( TEMP ( J ) * J  =  1  ♦  40 ) 

SUM=0.0 
DO  230  1=1*40 
SUM=SUM+TEMP( I ) 

230  CONTINUE 

SUM=SUM/40.0 
WRITE  ( LUW  *  240 )  SUM 

240  FORMAT  <//#»AVG#  TEMP  FOR  1  ST  SEC.  OF  RUN='*F9.3) 
WRITE  ( LUP  *  245 )  ( TEMP ( J )  ♦  J  =  1 *  40 ) 

245  FORMAT ( 10F8 • 3 ) 

CALL  CHAIN(BPLOT) 

END 
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****************  *****************  ****** 

*  * 

*  * 

*  PROCESS  PROGRAM  BPLOT  * 

*  - - — * 

*  * 

#  * 

*********  ***************************************  ****** 


DEFINE  FILE  5 1 ( 40 *80 *U » JO )  *  52 (40  *  80 *U ♦ JO ) 

DEFINE  FILE  53 ( 40  ♦  80 *U > JO ) 

DIMENSION  Z ( 40 ) »  X ( 40 ) 

COMMON  II  I  I (700) 

COMMON  LUR  *  LUW  *LUP*NDATA*NSS*NPUL*IPULW»DT 
COMMON  Y ( 400 ) 

I  DEV-7 

CALL  PLOTD(l) 

***  KEYBOARD  REQUEST  *** 

WRITE  ( LUW  #  10  ) 

FORMAT  (///♦ ‘PLOTTING  PROGRAM  *  TYPE-1  FOR  FLOW  PLOT*) 
WRITE  (LUW* 15) 

FORMAT ( 2 IX  * ' -2  FOR  TEMP  PLOT') 

CALL  FFINP(LUR*1 #0* IFILE* IEROR) 

I F I LE= I F I LE+50 

DETERMINE  RANGE  OF  ORDINATE  DATA  (Z) 

ZMAX=Q • 0 

ZM I N* 1000  •  0 

DO  40  JREC= 1  * 40 

READ  ( IFILE*  JREC  )  Z 

DO  35  1=1*40 

IF  (ZMAX-Z(I))  20*20*25 

ZMA X  =  Z (  I  ) 

GO  TO  35 

IF  (Z(I)-ZMIN)  30*30*35 
ZM I N=Z ( I ) 

CONTINUE 

CONTINUE 

WRITE  ( LUW  *45 )  ZMIN*ZMAX 

FORMAT  (5X»'DATA  VALUES  RANGE  BETWEEN '  * F 7  •  2  * ' ♦ *  •  F7 • 2 ) 
WRITE  ( LUW  *  50 ) 

FORMAT  ( 5X  * ' TYPE  ZMIN+ZMAX  DESIRED  FOR  PLOT ‘  ) 

CALL  FFINP  ( LUR *2  *  1 *ZM I N *  1 *ZMAX ♦ I EROR ) 
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CALL  APPROPRIATE  PLOTTING  SUBROUTINES 

DEFINE  SCALE  FACTORS  6  STARTING  PEN  LOCATION 

XSF=7. 2/1600.0 
ZSF=4.5/(ZMAX-ZMIN> 

XPPL*Q • 0 

ZPPL=ZMIN+(-3.0/ZSF) 

CALL  SCALF  ( XSF * ZSF * XPPL * ZPPL ) 

DEFINE  THE  PLOTTING  GRID 

***  X  AXIS  *** 

CALL  FGRID  ( 0 *0 . 0  * ZM I N * 100 . 0  *  16 ) 

***  Z  AXIS  *#* 

RLS (ZMAX-ZMIN) /10.0 

CALL  FGRID  ( 1 *0 . 0 *ZM I N * RL * 10 ) 

DEFINE  THE  DESIRED  ANNOTATION 

***  X  AXIS  *** 


XN=1600.0*0.09/7 .2 
ZN=ZMIN-(Q.25/ZSF) 

CALL  FCHAR (-XN*ZN* 0.09 *0.15 #0.0) 
WRITE  (  I  DEV ♦  60 ) 

60  FORMAT  ('  0 • * 18X * ' 10 ' *  1 8X * • 20  *  * 18X ♦ 
1 ' 30* *  18  X  * '40  •  ) 

XN=1600.0*3.0/7.2 
ZN  =  ZMIN-(0.5/ZSF  ) 

CALL  FCHAR (XN*ZN*0.15*0.2*0.0) 

WRITE  ( I  DEV ♦ 70 ) 

70  FORMAT  ('TIME  (SECS)*) 

C  ***  Z  AXIS  *** 

XN=1600. 0*0. 1/7.2 

CALL  FCHAR  (-XN*ZM IN *0.09 * 0.15  *1.57 ) 
WRITE  ( I  DEV  *  80 )  ZMIN.ZMAX 
80  FORMAT  ( F5 . 1  * 40X  •  F5 • 1 ) 
XN*1600.0*0.25/7.2 
ZN=ZMIN+( 1 .5/ZSF ) 

CALL  FCHAR (-XN.ZN# 0.15  *0.2.1.57 ) 

IF  (IFILE-51)  90*90*110 
90  WRITE  ( IDEV  *  100 ) 

100  FORMAT  ('FLOW  G.P.M.'J 
GO  TO  130 
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110  WRITE  ( I  DEV  *  120 ) 

120  FORMAT  ( » TEMP  DEG# F ' ) 

130  CONTINUE 

C  NOW  PLOT  THE  DATA 

CALL  FPLOT (+1*0#Q*ZMIN) 

CALL  FPLOT (+2*0#0*ZMIN) 

SUM=0 • 0 
DO  140  1=1*40 
140  X(  I  )=Q.Q 

DO  170  J= 1 *40 
DO  150  1=1*40 
X( I ) =SUM 
150  SUM=SUM+1 • 0 

READ  (IFILE'J)  Z 
DO  160  K=  1 *40 
XN=X ( K ) 

ZN  =  Z ( K  ) 

160  CALL  FPLOT ( 0  *  XN  *  ZN ) 

170  CONTINUE 

ZN=ZMIN-(3.0/ZSF) 

CALL  FPLOT  ( +1 *XN  *ZN ) 

C  ***  KEYBOARD  REQUEST  *** 

WRITE  ( LUW  *  180 ) 

180  FORMAT  ('TYPE-1  FOR  PLOT*-2  FOR  PUNCH*-3  FOR  EXIT') 
CALL  FFINP(LUR*1*0*IFLAG* IEROR) 

GO  TO  (190*200*260) *IFLAG 
190  WRITE  ( LUW  *19  5) 

195  FORMAT  ('TYPE-1  FOR  SCOPE*-2  FOR  PLOTTER') 

CALL  FFINP  (LUR* 1 *0* IJUMP  * IEROR ) 

GO  TO  ( 196*197) *  I  JUMP 

196  CALL  PLOTD(l) 

GO  TO  5 

197  CALL  PLOTD(O) 

GO  TO  5 

C  PUNCH  OUT  THE  DATA-USING  800  POINTS 

200  READ  (5*210) 

210  FORMAT  (10X) 

DO  230  J= 1 *40 
READ  ( 5 1 ' J )  X 

WRITE  (5*220)  ( X ( I  )  *  I  =  1  *  40  *  2 ) 

220  FORMAT  (10F7.3) 

230  CONTINUE 

DO  250  J= 1 *40 
READ  (52'J)  Z 
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WRITE  (5*240)  (  Z (  I  )  ♦  I  *  1 • 40  * 2 ) 

240  FORMAT  (10F8*3) 

250  CONTINUE 
260  CALL  VIAQ 
END 
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C  ******************************************** ********** 


c  *  * 

c  *  * 

C  *  PULSE  TESTING  ANALYSIS  PROGRAM  * 

C  *  * 

C  *  (PTAP)  * 

C  *  * 

C  *  * 

C  *  * 

C  *  TERRY  WILDMAN  * 


C  ********* ***** ****************************** ********** 

C 


5 

10 


COMMON  AREAP  * SSX  *SSY 

COMMON  I  TYPE » IOUT*WI •  WF * NNW » LUR * LUW * LUP 
COMMON  NN* ISS* IPUL*DTt I  QUAD* I  RED. I  EX  IT 
LUR=  1 
LUW=2 
LUP  =  6 

WRITE  ( LUW  *  5 ) 

FORMAT  ( ///>5X* • PULSE  TESTING  ANALYSIS  PROGRAM') 
WR I TE ( LUW  » 10 ) 

FORMAT  <  5X  * ' - - - «) 


C  ****KEYBOARD  REQUEST**** 


WRITE  <LUW»15) 

15  FORMAT ( // * 'ENTER- 1  TYPE* I  OUT* I  DEV* I  EX  IT ' ) 

CALL  FFINP  (LUR*4*0* ITYPE»0* IOUT*0* IDEV*0* IEXIT* IEROR) 

C  SELECT  READ  DEVICE 


IF  (IDEV-1)  20*20*40 


C  ****READ  DATA  VIA  KEYBOARD**** 


20  WRITE  ( LUW  *  25 ) 

25  FORMAT  ( // * ' ENTER-W I *WF *NNW • ) 

CALL  FFINP  (LUR*3*1*WI*1*WF*0*NNW*IEROR) 

GO  TO  (30*30*90*90) • I  TYPE 
30  WRITE  ( LUW  *  35 ) 

3  5  FORMAT  ( // * • ENTER-NN  *  I SS  *  I  PUL  * DT  *  I  QUAD  *  I  RED '  ) 

CALL  FFINP  (LUR*6*0*NN*0» ISS*0*IPUL*1*DT  *0*1  QUAD* 
/ 0  *  I  RED  *  I EROR ) 

GO  TO  95 

C  ****READ  DATA  VIA  CARD  READER**** 
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40  READ  (5*45)  WI*WF*NNW 
45  FORMAT  (2F10*4*I3) 

GO  TO  (50,50*90*90) *ITYPE 
50  READ (5*55)  NN *  I SS *  I  PUL »DT *  I  QUAD  *  I  RED 
55  FORMAT  (314, F7. 3*212) 

GO  TO  95 

90  CALL  LINK(WMAIN) 

95  CALL  L INK ( PMA IN ) 

END 
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***#*#*****-«•#■***#** ********* ************************** 
#  * 

*  * 

*  PROGRAM  WMAIN  * 

*  — - - . — --  * 

*  * 

*  * 

**##*#*#*  #-fc#*-tt*^**#*******#*#*-*tt********************** 


DIMENSION  W ( 80 ) *SWN(80) * AR ( 80 ) *  PH  I ( 8  0 ) 

COMMON  AREAP  * SSX  * SSV 

COMMON  I  TYPE  *  I  OUT »W I >WF #NNW*LUR #LUW*LUP 
COMMON  NN  » ISS ♦ I  PUL  » DT *  I  QUAD  *  I  RED  *  I  EX  I T 
DEFINE  FILE  51(5*80*U*ID) *52(5»8Q*U*ID) 

DEFINE  FILE  53 ( 5 . 80 *U *  I D ) *  54 ( 5  *  80 »U * ID ) 

C  ARE  FREQ.  RESULTS  ON  CARDS  OR  ALREADY  ON  FILE  $ 

IF  UTYPE-3)  10  #  1 0  •  30 

C  ON  CARDS 


10  DO  20  J= 1  *  NNW 

READ(5i  15)  W<  J) $ SWN ( J ) * AR ( J ) » PH  I ( J ) 

15  FORMAT  (4F10.4) 

20  CONTINUE 

C  NOW  WRITE  THE  RESULTS  IN  THE  APPROPRIATE  FILES 


30 


I  S  =  1 
I  E  =  40 
NEXT=1 

WRITE  ( 5 1 • NE  XT ) 
WRITE  ( 52 ' NEXT ) 
WRITE  ( 53 1  NEXT ) 
WRITE  ( 54  *  NE  XT ) 
CALL  LINK(APEN) 
END 


(  W  (  J  ) *  J= I S  > I E ) 

( SWN (J) * Js I S  » I E ) 
(AR< J) * Js I S • I E ) 
(PHI ( J) »J*IS»IE> 
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C  ****************************************************** 


C  *  * 

C  *  * 

C  *  PROGRAM  PMAIN  * 

C  *  * 

C  *  * 

C  *  * 


C  ****************************************************** 

c 

c 

c 

c 

c  the  following  options  are  provided 

c 

C  1  TIME  DATA  REDUCTION 

C 

C  2  QUADRATURE  METHOD 

C  A.  TRAPEZOIDAL 

C  B •  FILON'S 

C  C.  FAST  FOURIER  TRANSFORM 

C 
C 
C 

c 


COMMON  AREAP  *  SSX  * SSY 

COMMON  I  TYPE  * IOUT  *W  I *WF *NNW*LUR »LUW*LUP 
COMMON  NN#ISS»IPUL»DT*I QUAD  *IRED*IEXIT 
DEFINE  FILE  31 ( 5 *80 »U» ID ) *  52 ( 5 >80 *U* ID ) 
DEFINE  FILE  53 < 5  *  80  * U *  I D ) *  54 ( 5  *  80 *U *  I D ) 
DEFINE  FILE  5 5 ( 2 1  *  100  » U # I D ) t 56 ( 2 1 » 100  * U *  I D ) 
CALL  PREAD 
CALL  PWRIT 
IF  (IRED-1)  10*10*5 
5  CALL  REDN 
10  CONTINUE 
CALL  WVEC 
CALL  STEDY 

GO  TO  (15*20*25) *IQUAD 
15  CALL  L I NK ( TRAP ) 

20  CALL  LINMFILON) 

25  CALL  LINK(FFT) 

END 
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SUBROUTINE  PREAD 

DIMENSION  X( 1050 ) *Y( 1050) 

COMMON  AREAP  *  SSX  * SSY 

COMMON  I  TYPE » I OUT *WI »WF .NNW.LUR .LUW.LUP 
COMMON  NN.  I SS. I  PUL. DT. I  QUAD. I  RED. I  EX  IT 

C  ZERO  OUT  THE  VECTORS 

DO  5  J= 1 » 1050 
X< J)=0.0 
Y ( J ) =0 • 0 
5  CONTINUE 

C  IS  TIME  DATA  ON  CARDS  OR  ALREADY  IN  THE  FILES 

IF  { I  TYPE-1 )  10.10.40 
C  ON  CARDS 

10  READ  (5.15)  ( X ( J ) . J= 1  *  NN ) 

15  FORMAT  (10F7.3) 

READ  (5.20)  ( Y ( J  )  # J  =  1 *  NN ) 

20  FORMAT  (10F8.3) 


C  HOW  MANY  RECORDS  REQUIRED  $ 

NREC  = (NN/50 )+l 

C  WRITE  THE  TIME  DATA  ON  DISK 


IS=1 
I  E  =  50 

DO  35  NEX Ts 1 *NREC 
WRITE  ( 55 1  NEXT )  (X(J) .J-IS.IE) 

WRITE  ( 56 1  NEXT )  « Y< J) * J« IS. IE ) 

I SB I S  +  50 
IE- I E+50 
35  CONTINUE 
40  RETURN 
END 
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SUBROUTINE  PWRIT 
COMMON  AREAP • SSX  * SSY 

COMMON  I  TYPE • I OUT *W I  * WF * NNW * LUR * LUW * LUP 
COMMON  NN  *  I SS  *  I  PUL  *  DT *  I  QUAD  *  I  RED  *  I  EX  I T 
WRITE  ( LUW  *205 ) 

FORMAT  ( / / / *  1 8X  *  *  PULSE  ANALYSIS  PROGRAM’) 

WRITE  ( LUW  *210) 

FORMAT  (  1 8 X  *  *  * - - *  ) 

WRITE  ( LUW  *215)  NN 

FORMAT  ( /// *  5X  *  *  NO.  OF  DATA  POINTS  ='*I4) 

WRITE  ( LUW  » 220 )  ISS 

FORMAT  (  5  X  *  *  N  0  •  OF  S.S.  POINTS  = • *14) 

WRITE  ( LUW  *225)  DT 

FORMAT  ( 5X  * ' T  IME  INCREMENT  =  *  *  F7 • 3  *  ' SECS ’  ) 

WRITE  ( LUW  *230)  NNW 

FORMAT  C///*5X**NO.  OF  FREQ.  POINTS  CALCULATED- ’*13) 

GO  TO  (235*240*245) *  I  QUAD 
WRITE  ( LUW  *236) 

FORMAT  ( ///*5X* ' TRAPEZOIDAL  QUADRATURE  USED') 

GO  TO  250 
WRITE  ( LUW  *241 ) 

FORMAT  (///*5X*'FI LONS  METHOD  USED’) 

GO  TO  250 
WRITE  ( LUW  *  246 ) 

FORMAT  ( ///  *  5X  * • FAST  FOURIER  TRANSFORM*) 

CONTINUE 

IF  ( I  RED-1  )  280*280*275 
WRITE  ( LUW  *276 )  IRED 

FORMAT  (//*15X* ‘DATA  REDUCED-EVERY *  *  I  2  * • PO I  NT  CHOSEN*) 

CONTINUE 

RETURN 

END 
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SUBROUTINE  REDN 

DIMENSION  X ( 1050 )  * Y ( 1050 ) 

COMMON  AREAP.SSX *SSY 

COMMON  I  TYPE* I OUT  *WI * WF *NNW * LUR * LUW * LUP 
COMMON  NN* I SS* I  PUL* DT* I  QUAD. I  RED* I  EX  IT 

C  READ  RAW  TIME  DATA  INTO  CORE 

NREC= ( NN/50 ) +1 
I  S  =  1 
I  E  =  50 

DO  10  NEXT5*  1  * NREC 

READ  <55* NEXT )  ( X ( J ) * J= I S *  I E ) 

READ  <56* NEXT )  ( Y ( J ) * J= I S  *  I E ) 

I S= I S  +  5  0 
I E=  I E  +  50 
10  CONTINUE 

C  CALCULATE  THE  NEW  VECTOR  OF  DATA 

J=1 

DO  305  N= 1 *NN  *  I  RED 
X< J>=X<N> 

Y  (  J  )  *  Y  ( N  ) 

J  =  J+l 

305  CONTINUE 

C  CALCULATE  NEW  NUMBER  OF  DATA  POINTS 

NN  =  NN/ I  RED 

C  CALCULATE  NEW  NUMBER  OF  S.S.  POINTS 

I SSs I SS/ I  RED 

C  CALCULATE  NEW  NUMBER  OF  PULSE  POINTS 

I PUL= I PUL / I  RED 
C  CALCULATE  NEW  DT 

DT  =  DT *  I  RED 

C  WRITE  THE  REDUCED  DATA  ON  DISK 

NRECS (NN/50 ) +1 
IS=1 
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I  E  =  50 

DO  20  NEXT-1 *NREC 

WRITE  ( 55 • NEXT )  ( X ( J ) * j* I S *  I E ) 

WRITE  ( 56 1  NEXT )  ( Y ( J ) * J* I S • I E ) 

I Sa 1S+5Q 
I  E= I E  +  5  0 
CONTINUE 

WRITE  ( LUW  *310)  NN 

FORMAT  <//  *  5X  *  *  DATA  REDUCED  TO •  » 1 4  *  * PO I  NTS '  ) 
WRITE  ( LUW  *315)  ISS 

FORMAT  (  5X  *  1  NEW  NUMBER  OF  S.S.  POINTS  =  •  *  I  5 ) 
WRITE  ( LUW  *  320 )  IPUL 

FORMAT  { 5X  *  1  NEW  NUMBER  OF  PULSE  D I V  I S I ONS= *  *  I  5 ) 
WRITE  ( LUW »  32  5 )  DT 

FORMAT  ( 5  X  # 1  NEW  VALUE  OF  DT  =‘*F6.3) 

RETURN 

END 
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SUBROUTINE  WVEC 


CALCULATE  VECTOR  OF  FREQUENCY  POINTS 


DIMENSION  W ( 8 0 ) 

COMMON  AREAP iSSX  *SSY 

COMMON  I  TYPE  *  I  OUT  * WI * WF * NNW ♦ LUR $ LUW * LUP 
COMMON  NN* I SS ♦ I  PUL  *DT ♦ I  QUAD  » I  RED  * IEXIT 
WRITE  ( LUW  *401) 

401  FORMAT  (//*20X*»  FREQUENCY  RANGE') 

WRITE  (  L  U  W  *  4  0  2  )  WI*WF»NNW 

402  FORMAT  (5X*F10*4*'  TO  '*F10.5*‘  USING  '*16* 
1*  DIVISIONS') 

W I L=ALOG ( W I ) 

WFL=ALOG ( WF ) 

ND I V=NNW- 1 
DEL=(WFL-WIL) /NDIV 
DAL=DEL 
W( 1)*WI 

DO  405  Ms  2  *  NNW 
W(M)=EXP(DAL)*WI 
DAL=DAL+DEL 
405  CONTINUE 

WRITE  THE  FREQUENCY  VALUES  ON  DISK 


I  S  =  1 
I  E=40 
NEXT= 1 

WRITE  ( 5  1  ' NEXT )  ( W ( J ) * J= I S *  I E ) 

RETURN 

END 
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SUBROUTINE  STEDY 

DIMENSION  X ( 2 50 )  *Y(250) 

COMMON  ARE AP  *  SSX  * SSY 

COMMON  I  TYPE  * IOUT  *W I > WF » NNW * LUR * LUW * LUP 
COMMON  NN*ISS*IPUL*DT*I QUAD  *  I  RED  *  I  EX  IT 

C  READ  THE  S.S.  RAW  TIME  DATA  INTO  CORE 

NREC= (  ISS+IPUL) / 5  0+ 1 
I  S  =  1 
I  E  =  50 

DO  500  NEXT=1*NREC 

READ  (55*  NEXT )  ( X ( J ) * J= I S  t I E ) 

READ  (56'NEXT)  (  Y  (  J  )  *  J  =  I  S  .  I  E  ) 

I  S  = I S  +  50 
I  E=  I  E  +  50 
500  CONTINUE 

C  CALCULATE  THE  STEADY  STATES 

SUM=0 • 0 

DO  505  I=1*ISS 
SUM=SUM+X ( I ) 

505  CONTINUE 

SSX=SUM / I SS 

SUM=0 • 0 

DO  510  I  =  1  *  I  SS 
SUM=  SUM  +  Y ( I  ) 

510  CONTINUE 

SSY=SUM/ I SS 

WRITE  ( LUW  *515)  SSX 

515  FORMAT  (  /  / »  5  X  * •  I NPUT  S*S.  =**F9*4) 

WRITE  ( LUW  » 520 )  SSY 
520  FORMAT  (5X**OUTPUT  S.S.='fF9.4) 

DETERMINE  AREA  UNDER  PULSE  FOR  FREQ* 
CONTENT  CALCS 

AREAP=0 • 0 
I BGN= I SS+ 1 
I END= I BGN+ I  PUL 
DO  530  I  =  I  BGN  *  I  END 
AREAP=AREAP+(X( I )-SSX) 

530  CONTINUE 

AREAp  =  /\IRfEAP-(  (X(  I  BGN  )  +X  (  IEND)-2*0*SSX)/2*0) 
AREAP= ABS ( DT*AREAP ) 
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WRITE  ( LUW  *  540 )  AREAP 

FORMAT  (//  *  5X  *  *  AREA  UNDER  PULSE 

RETURN 

END 
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********* 

* 

* 

* 

* 

* 

* 

********* 


********************************************* 

* 

* 

PROGRAM  TRAP  * 

* 

* 

********************************************* 


DIMENSION  X ( 1050  )  *  Y ( 1050 ) 

DIMENSION  W(40) *AF(40) *BF<40) *CF(40) *DF(4Q) 
COMMON  AREAP  *SSX  *SSY 

COMMON  I  TYPE* I0UT*WI *WF *NNW *LUR *LUW*LUP 
COMMON  NN*ISS*IPUL*DT* I  QUAD* I  RED* I  EX  IT 
DEFINE  FILE  5 1 ( 5  ♦  80 *U *  I D ) *  32 ( 5  ♦  80 *U • I D ) 
DEFINE  FILE  53 ( 5 » 80 *U *  I D )  ♦  54 { 5  *  80 *U *  I D ) 
DEFINE  FILE  55 ( 2 1  *  100 *U *  I D ) *  56 ( 2 1  *  100 *U *  I D ) 

READ  THE  FREQ*  VECTOR  INTO  CORE 

NEXT  = 1 

READ  (51*  NEXT )  ( W < J ) * J= 1  *  40 ) 

READ  RAW  TIME  DATA  INTO  CORE 

NREC= ( NN/50 ) +1 
I S*  1 
I  E  =  50 

DO  600  NEXT= 1 »NREC 

READ  ( 55 1  NEXT )  ( X ( J ) * J* I S *  I E ) 

READ  C  5 6 • NEXT )  ( Y ( J ) * J* I S *  I E ) 

IS*IS+5Q 
I Es I E+5  Q 
600  CONTINUE 

CALC*  THE  QUADRATURE  PRODUCTS  FOR  ALL  THE 
FREQUENCY  POINTS 


DO  615  LL  =  1  *  NNW 

N.B.  -  QUADRATURE  CALCS  APPLY  ONLY  FROM  START 

OF  THE  PULSE 
I.E*  -  ISS+1 

A  =  0  •  0 
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B  =  0  •  0 
C  =  0.0 
D  =  0  •  0 

I  BGN  = I SS+ 1 
I END  = I BGN  + I PUL 

C  CALC  A  &  B 

K  =  0 

DO  605  I  =  I BGN  *NN 

A=A+ ( ABS ( Y ( I ) -SS  Y ) )*COS(W(LL)*K*DT) 

B  =  B+ ( ABS ( Y ( I )-SSY) ) *S I N ( W ( LL ) *K*DT ) 

K  =  K+1 

605  CONTINUE 
A= A*DT 
B=B*DT 

C  CALC  C  &  D 

K  =  0 

DO  610  I  =  I  BGN  *  I  END 

C=C+ ( ABS ( X ( I )-SSX) )*COS(W(LL)*K*DT) 

D*D+ ( ABS ( X ( I )-SSX) )*SIN(W(LL)*K*DT) 

K*K+1 

610  CONTINUE 

K  =  I  END- I BGN 

C  =  C-( ABS(X( I  END) -SSX) /2.0 ) *COS ( W ( LL ) *K*DT ) 
D=D- ( ABS ( X ( IEND) -SSX ) /2  •  0 ) *S I N ( W ( LL) *K*DT ) 
K  =  0 

C  =  C- ( ABS ( X ( I BGN) -SSX)/ 2.0 ) *COS ( W < LL ) *K*D T  ) 
D=D-( ABS(X( IBGN)-SSX)/2#0)*SIN(W(LL)*K*DT) 
C=C#DT 
D*D*DT 

C  STORE  THE  A*B*C*  AND  D  QUAD  PRODUCTS 

AF ( LL ) = A 
BF ( LL ) =B 
CF ( LL ) =C 
DF ( LL ) -D 
615  CONTINUE 

C  STORE  THE  RESULTS  IN  FILE  52 

NEXT= 1 

WRITE  ( 52 • NEXT )  ( AF ( J ) • J= 1 *40  ) 

NEXT=2 

WRITE  ( 52 1  NEXT )  ( BF ( J ) * J* 1 *40  ) 

NEXT- 3 

WRITE  (52 ' NEXT )  ( CF ( J ) •  Js 1 *40  ) 
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NEXT-4 

WRITE  ( 52 ' NEXT )  ( DF ( J ) * J- 1 *40  ) 

CALL  LINK(RSLTS) 

END 
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****************************************************** 


*  * 

*  * 

*  PROGRAM  FILON  * 

*  * 

*  * 

*  * 


****************************************************** 


DIMENSION  X( 1050  )  *Y(1050) 

DIMENSION  W(40) *AF(40) *BF(40) *CF(40) *DF(40) 
COMMON  AREAP  * SSX  * SSY 

COMMON  I  TYPE  *  I  OUT *W I  * WF  *  NNW  *  LUR  *  LUW  *  LUP 
COMMON  NN*ISS*IPUL*DT*I QUAD  *  I  RED  *  I  EX  IT 
DEFINE  FILE  5 1 ( 5 ♦ 80 *U *  I D ) • 5 2 ( 5 ♦ 80  * U . I D  ) 
DEFINE  FILE  53 ( 5 ♦ 80 *U *  I D  )  . 54 ( 5  *  80 *U *  I D ) 
DEFINE  FILE  55 ( 2 1  *  100 *U *  I D ) *  5 6 ( 2 1  *  100  * U *  I D ) 

READ  THE  FREQ.  VECTOR  INTO  CORE 

NEXT= 1 

READ  (51*  NEXT )  < W ( J ) * J= 1  *  40 ) 

READ  RAW  TIME  DATA  INTO  CORE 


NRECS ( NN/50 ) +1 
IS=1 
I  E  =  50 

DO  600  NEXT= 1 *NREC 

READ  ( 5 5  *  NEXT )  ( X < J ) * J* I S *  I E ) 

READ  (56*  NEXT )  < Y < J ) # J* I S • I E ) 

I S= I S+50 
IEME+50 
600  CONTINUE 

CALC  THE  QUAD  PRODUCTS  FOR  ALL  THE 
FREQUENCY  POINTS 


DO  650  LL  =  1  * NNW 

C  CHECK  SIZE  OF  THETA 

THETA=W (LL)*DT 
IF  (THETA-0.35)  605*605.610 
605  ALFA=(2.0*THETA**3/45.0)-(2.0*THETA**5/315.Q) 

ALFA= ALFA+ ( 2 . 0*THETA**7/472 5 • 0 ) 

BETA=(2.0/3.0)+(2.0*THETA**2/15.0)-(4.Q*THETA*#4 


if 
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3  •  B 


1/105.0) 

3ETA=BETA+(2.0*THETA**6/567.0 ) - ( 4 . 0*THET A**8 /2227 5 . 0 ) 
GAMMA® <4. 0/3.0 )-(2.0*THETA**2/15.0)+(THETA**4/270.0) 
GAMMA=GAMMA-( THETA**6/ 11340.0  >  +  ( THE TA**8/99 792 0.0 ) 

GO  TO  615 

610  ALFA= ( 1.0/THETA ) + ( S I N ( 2 . 0*THE TA ) / ( 2 . 0*THETA**2 ) ) 
ALFA=ALFA- ( (2.0* (SIN( THETA) ) **2 ) / ( THET  A**3 )  ) 

BETA= (  (COS (THETA)  ) **2  +  1 . 0 ) / ( THETA**2 ) 

BETA=2 • 0* ( BETA-SIN ( 2.0*THETA ) / ( THETA** 3 ) ) 

GAMMA® (4.0 /THETA**2)*( SIN (THETA) /THETA-COS (THET A) ) 

615  CONTINUE 

C  MAKE  NN  AN  ODD  NUMBER 

NN= ( (NN/2)*2)-l 
A  =  0  •  0 
B  =  0.0 
C  =  0.0 
D  =  0  •  0 
SI As0 • 0 
S2 A  =  0  •  0 
SI B  =  0  *0 
S2B=0 • 0 
S1D=0.0 
S2D®0  *0 
S1C=0*0 
S2C  =  0  *0 
I BGN® I SS+ 1 
JBGN= I SS+2 
I  END® I BGN+ 1 PUL  +  1 
JEND® IEND-1 
NNNN-NN-l 

CALC  A 
CALC  B 

K  =  0 

DO  620  I  ®  I BGN  *NN  *  2 

S1A®S1A+(ABS(Y( I ) -SSY ) ) *COS ( W ( LL ) *K*DT ) 
S1B=S1B+(ABS(Y< I ) -SSY ) ) *S I N ( W ( LL ) *K*DT ) 

K  =  K  +  2 

620  CONTINUE 
K=  1 

DO  625  J=JBGN*NNNN*2 

S2A=S2A+(ABS(Y(J) -SSY ) ) *COS ( W ( LL ) *K*DT ) 
S2B-S2B+(ABS(Y(J ) -SSY ) ) *S I N ( W ( LL ) *K*DT ) 

K  =  K+2 

625  CONTINUE 
C  CALC  A 
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K  =  0 

Z  =  ABS ( Y ( IBGN)-SSY) 

A= ALFA*Z*COS ( (WILL) *K*DT )+(3*l418/2*0) ) 

K  =  NN- I SS 
Z=ABS(Y<NN)-SSY) 

A=A-ALFA*Z-COS( ( W ( LL ) *K*DT )  +  (  3. 1418/2 . 0 ) ) 

A  =  A+BET  A*S1 A  +  GAMMA*S2A 
A=A*DT 

C  CALC  B 

K  =  0 

B=ALFA*ABS(Y( IBGN ) -SSY ) *COS ( W ( LL ) *K*DT ) 
K=NN-ISS 

B=B~ALFA*ABS ( Y { NN ) -SSY ) *COS ( W ( LL ) *K*DT ) 
B=B+BETA*S18+GAMMA*S2B 
B=B*DT 

CALC  C 
CALC  D 

K  =  0 

DO  640  1  =  1  BGN  *  I  END  *  2 

S1C=S 1C+ABS ( X ( I )-SSX)*COS(W(LL)*K*DT) 
S1D=S1D+ABS( X ( I )-SSX)*COS( W(LL)*K*DT) 

K  =  K+2 

640  CONTINUE 
K  =  0 

S1C=S1C-(ABS(X( I BGN ) -SSX ) ) *COS ( W ( LL ) *K*DT ) /2 . 
S1D=S1D-(ABS(X( I BGN ) -SSX ) ) *S I N ( W ( LL ) *K*DT ) /2 • 
K= I  END- I SS 

S1C=S1C-(ABS(X( I  END ) -SSX )  ) *COS ( W ( LL ) *K*DT ) /2 • 
S1D  =  S1D-(ABS(X( I  END ) -SSX )  ) *S I N { W ( LL ) *K*DT ) /2 . 

K  =  1 

DO  645  J= JBGN* JEND*2 

S2C=S2C+( ABS(X( J ) -SSX ) ) *COS ( W ( LL ) *K*DT ) 
S2D=S2D+(ABS(X( J)-SSX) ) *S I N ( W ( LL ) *K*DT ) 

K  =  K+2 

645  CONTINUE 
C  CALC  C 

K  =  0 

Z= ABS ( X ( IBGN )-SSX ) 

C=ALFA*Z*COS ( (WILL) *K*DT )+(3.l4l8/2*0) ) 

K=  I  END- I SS 
Z  =  ABS ( X ( I  END ) -SSX ) 

C=C-ALFA*Z*COS ( (W(LL)*K*DT)+(3. 1418/2*0) ) 
C=C+BETA#S1C+GAMMA*S2C 
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C=C*DT 


K  =  0 


C  CALC  D 


K  =  0 

D= ALFA*ABS ( X ( IBGN ) -SSX ) *COS ( W ( LL ) *K*DT ) 

K=  IEND-ISS 

D=D-ALFA#A8S ( X ( I  END) -SSX ) *COS ( W  (  LL ) *K*DT ) 

D=D+BETA*S1D+GAMMA*S2D 

D=D*DT 

C  STORE  THE  A*B*C*  AND  Dt  QUAD  PRODUCTS 

AF ( LL ) =A 
BF ( LL ) =B 
CF ( LL ) s C 
DF ( LL ) =D 
650  CONTINUE 


C  STORE  THE  RESULTS  IN  FILE  52 


NEXT= 1 

WRITE  (52‘NEXT) 
NEXT  =  2 

WRITE  ( 52  *  NEXT ) 
NEXT=  3 

WRITE  (52'NEXT) 
NEXT= A 

WRITE  (52‘NEXT) 
CALL  LINK(RSLTS) 
END 


( AF ( J ) » J= 1*40) 
( BF ( J ) *  J= 1 » 40  ) 
(CF< J) *  J= 1  *  40  ) 
( DF ( J ) *  J* 1 # 40  ) 


Ti *3  =  0 
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****************************************** ************ 


*  * 

*  * 

*  PROGRAM  FFT  * 

*  - - ---  * 

*  * 

*  * 


****************************************************** 


DIMENSION  I  SUB ( 40 ) 

DIMENSION  AR( 1050) *AI ( 1050) 

DIMENSION  W(40) »AF(40) *BF(40) »CF(40) *DF(40) 
COMMON  AREAP * SSX  *SSY 

COMMON  ITYPE*IOUT  *  W  I *  WF *NNW*LUR*LUW*LUP 
COMMON  NN* ISS* I  PUL t DT *  I  QUAD  *  I  RED# I  EX  IT 
DEFINE  FILE  5 1 ( 5  *  80 *U *  I D  )  *  52 ( 5  *  80 *U *  I D  ) 
DEFINE  FILE  53 ( 5  *80 *U* ID >  *54 ( 5  *  80 *U* ID ) 
DEFINE  FILE  55 < 2 1  *  100 *U *  I D ) *  56 ( 2 1  *  100 *U » I D ) 


SPECIFY  THE  PARAMETERS  FOR  TNFRM 


N  =  NN 
I  N  V-  1 
I  NP=  1 

CALCULATE  NM 


NPROD-N 
NM  =  Q 

5  NM=NM+1 

NPRODsNPROD/2 
IF  (NPROD-1)  10*10*5 
10  CONTINUE 

CALCULATE  A  NEW  FREQUENCY  VECTOR 

SO  THAT  THE  RESULTS  ARE  SUITABLE  FOR  PLOTTING 


K  =  0 
1  =  1 
L=  1 

DO  20  J= 1  * NNW 
W( J)=6.283*J/(NN*DT> 
I  SUB ( J ) - J 
20  CONTINUE 


C 


WRITE  THE  NEW  FREQ  VECTOR  BACK  ON  THE  DISK 
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NEXT= 1 

WRITE  ( 5 1  *  NEXT )  ( W ( J ) * j= 1  * 40 ) 

C  READ  THE  INPUT  TIME  DATA  INTO  CORE 

NREC* ( NN/ 50 ) +1 
15  =  1 
I  E  =  50 

DO  30  NEXT=1*NREC 

READ  (55*  NEXT )  < A  I ( J ) * J« I S * I E ) 

15=  I S+5  0 
I  E= I E  +  5  0 
30  CONTINUE 

SUBTRACT  THE  S»S.  VALUE  FROM  THE  TIME  DATA 
AND  STORE  IN  AR ( K ) 

DO  35  J* 1  *  NN 
AR  (  J ) = ABS ( A  I ( J)-SSXJ 
35  CONTINUE 

CALL  THE  FAST  FOURIER  TRANSFORM 

CALL  TNFRM(N*NM*AR*AI * INV* INP ) 

SELECT  THE  APPROPRIATE  A  AND  B  VALUES  FROM 
THE  VECTORS  AR*  AND  A I  *  RESULT  I NG  FROM  TNFRM 

DO  40  J= 1  * NNW 
I  =  I  SUB ( J  ) 

CF ( J ) = AR  (  I  ) 

DF ( J ) s A  I  (  I  ) 

40  CONTINUE 

READ  THE  RESPONSE  DATA  INTO  CORE 

I  S  =  1 
I  E  =  50 

DO  45  NEXT=1*NREC 
READ  ( 56  *  NEXT )  ( A  I ( J ) * Js I S *  I E ) 

I  S= I S  +  50 
I  Es  IE  +  50 
45  CONTINUE 

SUBTRACT  THE  S«S#  VALUE  FROM  THE  TIME  DATA 
AND  STORE  IN  AR ( K ) 

NN  =  N 

DO  50  J= 1  *  NN 
AR ( J ) =ABS ( A  I ( J)-SSY) 
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50  CONTINUE 

C  CALL  THE  FAST  FOURIER  TRANSFORM 


CALL  TNFRM(N*NM*AR*AI » INV» INP ) 

C  SELECT  THE  APPROPRIATE  C  AND  D  VALUES 

DO  55  J  =  1  * NNW 
I  =  I  SUB ( J ) 

AF ( J ) = AR (  I  ) 

BF ( J ) aA I (  I  ) 

55  CONTINUE 


C  STORE  THE  RESULTS  IN  FILE  52 


NEXT=  1 

WRITE  ( 52 ' NEXT ) 
NEXT=2 

WRITE  ( 52 1  NEXT ) 
NEXT  =  3 

WRITE  ( 52 ' NEXT ) 
NEX  T  =  4 

WRITE  ( 52 1  NEXT ) 
CALL  LINK(RSLTS) 
END 


( AF ( J ) ♦ J= 1 *40  ) 
( BF ( J) » J= 1 *40 ) 
( CF ( J ) *  J* 1  *  40 ) 
( DF ( J ) *  J= 1 *40  ) 
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£******  ***************  *********  *  *  *  *  *  *  ********  *  ************* 
C  **  SUBROUTINE  TNFRM  ** 

£** ***************************************** *************** 

C 

C  SUBROUTINE  ‘TNFRM'  CALCULATES  THE  FOURIER  TRANSFORM 

C  USING  THE  COOLEY-TUKEY  TRANSFORM  AS  WRITTEN  BY 

C  TRINITY  UNIVERSITY  COMPUTER  CENTER. 

C 

c 

c********* ************************************ ************* 
c  **  PARAMETERS  ** 

C 

C  N  THE  NUMBER  OF  DATA  POINTS  WHICH  MUST  BE  A  POWER 

C  OF  2. 

C*****NOTE.  MAXIMUM  DIMENSION  IS  TO  BE  1024  UNLESS  THE 


DIMENSION  OF  ACOS  IS  CHANGED  IN  THE  PROGRAM. 

NM  -  THE  EXPONENT  OF  2t  I.E.  N=2**NM 

AR  -  THE  REAL  INPUT  AND  OUTPUT  ARRAY »  WHICH  HAS  THE 
DIMENSIONS  OF  N 

A I  -  WHICH  MAY  BE  THE  IMAGINARY  INPUT  BUT  IS  ALWAYS 
THE  IMAGINARY  OUTPUT  ARRAY 

INV  -  IF  YOU  GET  THE  TRANSFORM »  IF  =2*  YOU  GET  THE 

INVERSE  TRANSFORM 

INP  -  IF  =1*  IT  INDICATES  REAL  INPUT  DATA  ONLY*  IF  =2* 
IT  INDICATES  COMPLEX  INPUT  ONLY 


C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 

£*************■&•*****#■###  *#■&***#*  *************************** 
c  **  GENERAL  COMMENTS  ** 

C 

C  THIS  LISTING  FOR  'TNFRM'  HAS  BEEN  USED  AS  RECIEVED 

C  WITHOUT  MODIFICATION  BY  THE  THE  AUTHOR  OF  THIS 

C  THESIS. 

C 

c 

^*****^************************************-i*,************-H-** 

c 

c 


SUBROUTINE  TNFRM  ( N *NM * AR *  A I  *  I NV * INP ) 

C*****NOTE.  ACOS  MUST  ALWAYS  BE  DIMENSIONED  TO  NX4 • 

DIMENSION  ACOS (257) 

DIMENSION  AR ( 1 )  *  A  I  ( 1 ) 

COMMON  ARE AP  *  SSX  *  SSY 

COMMON  ITYPE* I OUT  »WI *WF»NNW*LUR»LUW*LUP 
COMMON  NN*ISS*IPUL*DT*I QUAD  *  I  RED  *  I  EX  I T 
NX  1  =  N /2 
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NX2=NX 1 + 1 
NX  3  =  N /4 
NX4=NX3+1 
PI  2  =  6.2831833072 
DO  3  1  =  1* NX4 
I  P=  I - 1 

ACOS (  I  )=COS  ( P I  2* I P /N ) 

3  CONTINUE 

PERFORM  FIRST  TWO  POINT  TRANSFORM. 
GO  TO  40  IF  INPUT  DATA  COMPLEX. 

I  F (  I NP-2 )  30  *  40*  30 
30  DO  4  J=1*NX1*1 
NXO= J+NX1 
ARQ  =  AR ( J ) 

AR(J)=AR(J)+AR(NXC) 

AR ( NXO ) =ARQ-AR ( NXO ) 

4  CONTINUE 
GO  TO  56 

40  DO  50  J  =  1#  NX  1 
NXO= J  +  NX  1 
ARQ= AR ( J ) 

AR ( J) = AR ( J ) +AR ( NXO ) 

AR ( NXO ) = ARQ-AR ( NXO) 

A  I  Q  =  A  I  (  J  ) 

A  I  ( J )  =  A  I  ( J )  +  A  I (NXO) 

A  I (NXO)  =  A  I Q  -  A  I  (NXO) 

50  CONTINUE 
56  CONTINUE 

DO  6  L  =  2  »  NM  *  1 
NL=NM-L 
N X  =  2  *  * ( NL ) 

NZ=NX+NX 
NN=N-NZ+1 
DO  7  K= 1  *  NN  *  NZ 
I S=K“ 1 
IX= IS/NX 

C  CALL  SUBROUTINE  FOR  BIT  REVERSAL 

CALL  REVER  (NM*IX.IXR) 

IR=  IXR+1 

I F (  IXR-NX3 ) 9  *9  *  8 

8  JH= IR-NX3 
JG=NX2-IXR 
C=“ACOS ( JG ) 

D= ACOS ( JH ) 

GO  TO  10 

9  JF=NX4-IXR 
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C= ACOS ( IR ) 

D= ACOS ( JF ) 

GO  TO  10 

10  I F (  I  NV-2 ) 23 ♦  22*  23 

22  D  =  -D 

23  CONTINUE 

DO  11  M=  1  *  NX*  1 
NJ=  IS+M 
NX J-N J+NX 
A  I Q= A  I  (NJ  ) 

ARQ=AR ( N J ) 

I  F  (  I  N  P-2  )  63,  .13*  65 

65  I F ( L-2  )  13*  12*  13 

12  CONTINUE 

AR(NJ)  =  AR ( NX J ) *C  +  AR(NJ) 

A  I  ( N J )  =  AR ( NX J ) *D 

AR ( NXJ ) =-AR ( NJ )  +  ARQ  +  ARQ 

A  I  ( NX  J ) =- A  I  (NJ) 

GO  TO  14 

13  AR(NJ)  -  AR ( NX J ) *C  -  AI(NXJ)*D  +  AR ( N J ) 
AI  (NJ)=AR(.NXJ)*D  +  A  I ( NX J ) *C  +  AI(NJ) 

AR ( NX J ) =- AR ( NJ ) + ARQ+ARQ 

A I  (NXJ)  =  -AI(NJ)  4  A  IQ  +  A I Q 

14  CONTINUE 

11  CONTINUE 
7  CONTINUE 
6  CONTINUE 

DO  1  I  =  1  *  N 
I X= I - 1 

CALL  SUBROUTINE  FOR  BIT  REVERSAL 

CALL  REVER  ( NM *  I  X • I XR ) 

IXXR= IXR+1 
IF (  IXXR-I  ) 1 » 1 *101 
101  CONTINUE 

ARQ= AR ( IXXR ) 

AR (  IXXR  )  = AR ( I ) 

AR (  I  ) = ARQ 
A  I Q= A  I  (  IXXR) 

AI  (  IXXR  )  = A  I  ( I  ) 

A  I  (  I  )  =  A  I Q 
1  CONTINUE 

I  F (  I  NV-2 ) 70*66*70 

66  PD=1.0/N 

DO  70  I =1 *N 
AR (  I  ) = AR (  I )*PD 
A  I  (  I  )  = A  I  (  I  )*PD 
70  CONTINUE 
RETURN 
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* 

#****- *******  *****#*##■**##****■*•* 

1  *********  *  * 

*  SUBROUTINE  REVER  * 
*#**************#*********-***■#•#****#*#*■* 

1*********** 

* 

*  'REVER*  CARRIES  OUT  THE  BIT  REVERSAL 

1  OPERATION 

*  REQUIRED  BY  SUBROUTINE  ' TNFRM ' • 

*  IT  WAS  USED  AS  SUPPLIED  BY  THE  TRINITY 

*  UNIVERSITY  COMPUTER  CENTER. 

* 

*  *  ***********  *  *  *  *  *  *  *  *******  ******  *  *****  * 


1*********** 

* 

*  CALL  REVER  (NM*IX.IXR) 

* 

*  NM  -  POWER  OF  TWO  (NUMBER  OF  BITS) 

*  IX  -  NUMBER  TO  BE  REVERSED 

*  IXR  -  NUMBER  REVERSED 

* 


* 

.a.****.#**.********************#*********** 


I*********** 

* 

* 


ENT 

REVER 

REVER 

DC 

0 

LDX 

I  1 

REVER 

LD 

I  1 

+ 1 

LOAD  IX 

STO 

NUMBR 

STORE  AS  NUMBR 

LD 

11 

0 

LOAD  NM 

STO 

REVRS 

LDX 

I  1 

REVRS 

INTO  XR1 

SLT 

32 

CLEAR  A  AND  Q 

1 

REGISTERS 

STO 

REVRS 

SHIFT 

LD 

NUMBR 

SRT 

1 

SHIFT  A  AND  Q  RIGHT 

1 

ONE 

STO 

NUMBR 

LD 

REVRS 

LOAD  REVERSE  SUM 

1 

INTO  A 

SLT 

1 

SHIFT  Q  LEFT  INTO  A 

1 

BY  ONE 

STO 

REVRS 

MDX 

1 

-1 

DECREMENT  XR1 

MDX 

SHIFT 

* 
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Q  ********  *  ***********************  *******  ****  *********** 


c  *  * 

C  *  * 

C  *  PROGRAM  RSLTS  * 

C  *  * 

c  *  * 

Q  *  * 


c  ****************************************************** 

c 

c 

DIMENSION  AF ( 40 )  *BF(40)  #  CF ( 40  )  »DF(4Q) 

DIMENSION  VH  80)  *SWN(80 )  * AR ( 8  0 ) *  PH  I ( 8  0 ) 

COMMON  AREAP  * SSX  »  SSY 

COMMON  I  TYPE  *  I  OUT  tWl #WF*NNW*LUR  *LUW*LUP 
COMMON  NN  » I SS ♦ I  PUL  *  DT  #  I  QUAD* I  RED# I  EX  IT 
DEFINE  FILE  5  1  ( 5  *  80 *U *  I D )  ♦  52 ( 5 » 80 *U *  I D ) 

DEFINE  FILE  5 3 ( 3  ♦  80  * U *  I D )  » 54 ( 3 # 80 » U » I D  ) 

C  READ  THE  FREQUENCY  VECTOR  INTO  CORE 

NEXT=  1 

READ  ( 5 1 1  NEXT )  ( W l J ) ♦ J= 1 » 40 ) 

READ  THE  QUADRATURE  PRODUCTS  OFF  FILE  52 
AND  INTO  CORE 

NEXT=  1 

READ  (52*  NEXT )  (AF(J) #J=1*40) 

NEXT  =  2 

READ  (  52  ' NEXT )  ( BF ( J ) # J= 1 *40 ) 

NEX  T  =  3 

READ  (52'NEXT)  ( CF ( J ) # J= 1 *40 ) 

N  E  X  T  =  4 

READ  (52'NEXT)  ( DF ( J )  * J=  1 » 40 ) 

CALC  THE  AMP  RATIO  PHASE  LAG  AND  FREQUENCY 
C  CONTENT  FOR  ALL  FREQUENCY  POINTS 


DO  73  5  L L  =  1  *  NNW 
A=AF ( LL ) 

B=BF ( LL ) 

C  =  CF ( LL  ) 

D=DF( LL ) 

C  CALC  REAL  PART 
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RE=(A*C+B*D)/ (C**2+D**2 ) 

C  CALC  IMAG.  PART 

AIM=  (  A*D-B*C  )  / (C**2+D**2  ) 

C  CALC  AMP  RATIO 

AR ( LL )=SQRT ( RE**2+A I M**2 ) 

C  CALC  PHI 

IF  (RE)  710.710*705 

705  PHI <LL)=ATAN( AIM/RE) *180. 0/3. 14 18 
GO  TO  730 

710  IF  (AIM)  715.715*720 

715  PHI (LL)=( ATAN(ABS( AIM/RE) )*180.0/3.1418)-180.0 
GO  TO  730 

72  0  PHI  ( LL  )=- 180.0-A TAN (ABS( AIM/RE)  )*18Q. 0/3. 1418 
730  CONTINUE 

C  CALC  FREQ  CONTENT 

SWN(LL)=SQRT(C**2+D**2 ) /AREAP 
735  CONTINUE 

C  WRITE  THE  FREQUENCY  VECTORS  ON  DISK 


15=1 
IE  =  40 
NEXT=  1 

WRITE  ( 5 1 1  NEXT ) 
WRITE  (52'NEXT) 
WRITE  ( 53 • NEXT ) 
WRITE  ( 54 1  NEXT ) 
CALL  LINK(APEN) 
END 


(W(j)  » J=  I  S  ♦  I E ) 

( SWN (J) *J=IS*IE) 
( AR( J) *J=IS» IE) 
(PHI (J) . J= I S  *  I E ) 
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*  MAINLINE  APEN  * 
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***********  * ****** *  * *****************  *  *  *  ************** 


DIMENSION  W< 80 ) * SWN ( 80 )  * AR ( 80 ) ♦ PH  I  ( 80 ) 
COMMON  AR  EAP  »  SSX  *  SSY 

COMMON  I  TYPE  *  I  OUT »W I  * WF * NNW * LUR ♦ LUW » LUP 
COMMON  NN *  I SS  *  I  PUL  *  DT  *  I QUAD* I  RED  » I  EX  I T 
DEFINE  FILE  5 1 ( 5  ♦  80 *U ♦ I D )  *  5 2 { 5  *  80 *U *  I D ) 
DEFINE  FILE  5 3 ( 5  *  80 »U ♦ I D )  *  54 ( 3  *  30  * U *  I D ) 
GO  TO  ( 800  *825 *800 )  ♦  IOUT 

READ  THE  FREQUENCY  VECTORS  INTO  CORE 


800 


805 


810 


815 

820 

825 

830 

835 

840 


CONTINUE 
I  S  =  1 
I  E  =  40 


NEXT=1 

READ  (51'NEXT) 
READ  (52‘NEXT) 
READ  (53*  NEXT ) 
READ  ( 54 • NEXT ) 
WRITE  ( LUP  *305 ) 
FORMAT ( ///  ♦  10X  * 


(  W  (  J  ) *J=IS*IE) 

( SWN ( J )  * J= I S  *  I E  ) 

( AR ( J )  * J=  IS  *  IE ) 

( PHI  ( J ) * J= IS.  IE ) 

•FREQUENCY •  *  5X ♦ ' FREQUENCY 


1 •AMPLITUDE' *  5X  *• PHASE '  ) 


*  5  X  * 


WRITE  ( LUP  *  8 1 0 ) 

FORMAT  ( 1 0  X  ♦  '  RAD/SEC  •*5X**  CONTENT  *  *  5  X  *  '  RATIO 
15X  * • ANGLE •  ) 

DO  820  L  = 1  ♦  N N W 

WRITE  (LUP *8 15)  W ( L ) * SWN ( L ) * AR ( L )  » PH  I ( L ) 

FORMAT  (  10X*F9.4*5X*F9«4*3X*F8.3*5X*F9*3) 

CONTINUE 

IF  ( I  OUT  —  2  >  830  *  825  *825 
CALL  LINK(APLOT) 

IF  (IEXIT-1)  835*835*040 
CALL  LINK(PTAP) 

CALL  EXIT 
END 
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MAINLINE  APLOT  * 

- — * 

* 

* 

*  ************************  ************************ 


DIMENSION  X(40) *Y(40) *LABX( 10 ) *LABY(10) 
COMMON  AREAP  *  SSX  *  SSY 

COMMON  I  TYPE* IOUTtWI ♦ WF * NNW ♦ LUR * LUW * LUP 
COMMON  NN* I SS* I  PUL* DT *  I  QUAD* I  RED* I  EX  IT 
DEFINE  FILE  5 1 ( 5  *  80  * U *  I D ) *  52 ( 5  *  80  * U *  I D ) 
DEFINE  FILE  5 3 { 5  ♦  80 *U *  I D )  *  54 ( 5 . 80 ♦ U *  I D ) 


C  PREDETERMINED  PARAMETERS 


N  =  40 
I  S  =  2 
I  G  =  2 
I  T  =  1 
13  =  1 
I  R  =  2 
XM I N  =  W  I 
XMAX= WF 
NF I LX  =  5  1 
NPT  S=NNW 


C  CALC  NO  OF  LOG  CYCLES  ON  X  AXIS 

5  PROD=XM I N#1 • 0 1 

NC YCX=0 

10  NCYCX=NCYCX+1 
PROD= 10 .0*PROD 
IF  (XMAX-PROD)  20*20*10 
20  CONTINUE 

WRITE  ( LUW  *  3  0 ) 

30  FORMAT  ( 10X * ' PLOTT ING  PARAMETERS') 

WRITE  ( LUW  *  40 ) 

40  FORMAT  (  ' ENTER- I D *  I L ♦ I U *  I P * NF I LY * YM I N ♦ YMAX '  ) 

CALL  FFINP  ( LUR *  7 *0  *  I D *0  *  I L *0  * IU *0 ♦ I P *  0 *NF I LY *  1  * 
1YMIN* 1 » YMAX* IEROR ) 
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C  IF  Y  LINEAR  DO  NOT  CALC  NO  OF  LOG  CYCLES 

IF  (IL-1)  45*60*45 
45  IF  (IL-3)  48*60*48 

C  CALC  NO  OF  LOG  CYCLES  ON  Y  AXIS 

48  PROD=YMIN*1#01 
NC YC Y=Q 

50  NCYCY=NCYCY+1 
PROD=iO.O*PROD 
IF  (YMAX-PROD)  60*60*50 
60  CONTINUE 

C  SPECIFY  THE  X  LABEL 

WRITE  ( LUW  *  70 ) 

70  FORMAT  ('ENTER  THE  LABEL  FOR  THE  X  AXIS') 

CALL  FFINP  ( LUR *  1  *  102 *LABX *  I EROR ) 

C  SPECIFY  THE  Y  LABEL 

WRITE  (LUW *80) 

80  FORMAT  ('ENTER  THE  LABEL  FOR  THE  Y  AXIS') 

CALL  FFINP  (LUR*1 . 102.LABY .IERQR) 

C  CALL  THE  PLOTTING  SUBROUTINE 

CALL  NWPLT  ( X  * Y * N  ♦  I D ♦ I S *  I G » I L *  I T *  I B *  I U *  I R ♦ I P * XMAX * 

1 XM I N  »  NC YCX *  YMAX *  YM I N * NC YC Y  * LABX  *  LAB Y * NF I  LX  * NF I LY  * NPTS ) 

C  CHECK  WHETHER  LINK  OR  EXIT 

IF  (IEXIT-1)  90*90*100 
90  CALL  LINK(PTAP) 

100  CALL  EXIT 
END 
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